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EXECUTIVE  SUMMARY 


This  final  report  presents  research  carried  out  by  Berkeley  Research  Associates,  Inc. 
under  Contract  N00014-85-C-0680  with  the  Naval  Ocean  Research  and  Development  Activ¬ 
ity  (NORDA).  The  report  covers  the  performance  period  3  July  1985  through  30  September 
1987.  During  this  period,  in  conjunction  with  our  ongoing  effort  with  NORDA,  we  carried 
^out  research  which  would  improve  the  accuracy  of  the  NPE  (Nonlinear  Progressive- Wave 
Equation),  TOPS  (Thermodynamic  Ocean  Prediction  System)  and  Hibler  Sea  Ice  models 
and  would  reduce  the  discrepancies  between  the  models  and  actual  experimental  observa¬ 
tions.  In  the  following  paragraphs  we  will  summarize  the  work  carried  out  on  these  models 
during  this  performance  period. 

In  order  to  simulate  a  greater  variety  of  ocean  environments,  the  NPE  computer  code 
has  been  modified  to  include  a  damping  region.  The  inclusion  of  these  enhancements 

r  - 

has  provided  a  means  to  better  understanding^real  oceans.  Simulating  only  the  region 
surrounding  and  encompassing  a  moving  pulse  by  moving  in  the  frame  of  reference  of  that 
pulse,  the  NPE  contains  a  boundary  at  the  bottom  of  its  simulation  grid  capable  only 
of  reflecting  the  propagating  sound  waves.  In  a  real  ocean  environment,  one  finds  that 
sound  waves  propagate  downward  in  depth  and  reflect  from  such  boundaries  as  the  ocean 
sediment,  yet  they  do  not  completely  reflect  from  these  boundaries  as  they  would  from  a 
hard  boundary.  Rather,  some  of  the  sound  waves  will  penetrate  the  bottom,  reflecting, 
refracting,  diffracting  and  steepening  within  the  sediment,  some  upward  toward  the  surface, 
some  downward  into  the  sediment  (and  out  of  the  simulation  region).  Therefore,  a  damping 
region  is  required  to  prevent  sound  waves  from  reflecting  from  the  bottom  of  the  grid  and 
reentering  the  simulation.  We  have  therefore  included  a  damping  region  into  the  NPE 
code  and  have  been  successful  in  eliminating  these  undesirable  relections. 

Also  important  in  allowing  for  the  simulation  of  a  greater  variety  of  ocean  environ¬ 
ments  is  the  implementation  of  a  range  dependent  sound  speed  profile.  As  the  NPE 
demarcates  the  boundary  between  an  ocean  and  its  sediment  by  a  jump  in  sound  speed, 
the  ocean  with  one  sound  speed  and  the  sediment  with  a  somewhat  higher  one,  it  is  easy 


to  see  that  a  sound  speed  profile  that  is  constant  in  range  will  allow  for  only  flat  bound¬ 
aries  between  the  ocean  and  the  sediment.  By  incorporating  a  sound  speed  profile  which 
varies  with  range,  more  complicated  ocean-sediment  boundaries  may  be  simulated.  We 
have  therefore  included  a  function  into  the  NPE  code  which  is  called  each  timestep  to 
change  the  ocean-sediment  boundary  based  on  the  range  of  the  simulation  grid  and  have 
been  successful  in  simulating  various  ocean  systems  with  sloping  bottoms. 

While  it  is  important  to  increase  the  variety  of  problems  being  simulated  by  the 
NPE,  it  is  also  important  that  we  continue  to  improve  upon  its  accuracy.  For  some  time 
now,  the  NPE  code  has  maintained  a  timestep  constraint  which  adjustes  the  size  of  the 
timestep,  based  on  the  nonlinearity  in  the  simulation,  to  an  accuracy  dictated  by  the  user. 
Recently,  the  NPE  code  has  been  modified  to  also  include  a  diffusion  timestep  constraint. 
By  writing  the  diffraction  portion  of  the  NPE  in  the  form  of  a  continuity  equation,  the 
root  mean  square  velocity  is  calculated  and  used  to  obey  the  Courant-Friedrichs-Lewy 
condition  within  a  specified  level  of  accuracy.  The  use  of  this  constraint  is  an  improvement 
over  the  previous  method  used,  namely  that  of  manually  calculating  the  initial  timestep 
and  keeping  it  constant  throughout  the  (linear  portion  of  the)  simulation. 

Also  used  to  provide  for  greater  accuracy  is  the  implementation  of  a  tensile  cutoff. 
In  the  “real”  world,  one  finds  that  pressures  that  are  less  than  zero  do  not  exist.  In  a 
“simulated”  world,  one  must  “tell”  the  simulation  that  such  pressures  do  not  exist.  We 
have  therefore  included  a  tensile  cutoff  into  the  NPE  code  which  prevents  the  pressure 
from  falling  below  zero.  By  using  an  assumed  adiabatic  equation  of  state,  we  calculate 
what  the  overdensity,  pmin  must  be  to  correspond  to  a  zero  total  pressure.  We  then  solve 
the  NPE  as  usual  and,  for  each  timestep,  set  any  overdensity  falling  below  pm,n  to  the 
value  of  pmin ■  Therefore,  the  total  pressure  at  any  point  or  time  of  the  simulation  will  be 
at  least  zero. 

Of  course  the  NPE  requires  that  an  initial  pulse  with  certain  initial  conditions  be 
input  into  the  code.  While  in  the  past  we  have  primarily  used  simple  initial  pulses  for 
diagnostic  purposes  (being  that  it  is  easier  to  predict  what  a  simple  pulse  will  do  as  it 
propagates  than  to  do  so  for  a  complicated  pulse),  it  is  important  for  the  NPE  to  have 
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some  way  of  including  a  real  ocean  initial  condition.  We  have  therefore  written  the  code 
required  to  take  the  output  of  the  nonlinear  CSQ  code  from  Sandia  and  input  it  into  the 
NPE  as  an  initial  condition.  The  NPE  will  therefore  be  capable  of  simulating  oceans  in 
which  the  initial  conditions  are  indicative  of  a  real  ocean. 

The  remainder  of  the  NPE  section  of  this  report  is  used  to  show  some  results  which 
have  been  accumulated  throughout  this  performance  period.  Shallow,  deep  and  ice  con¬ 
taining  oceans  are  simulated,  and  the  results  of  error  and  diagnostic  tests  are  shown. 

In  predicting  temperature  profiles  in  the  upper  ocean,  different  statistical  models  may 
be  used.  Of  those  presently  used,  the  simplest  is  that  of  climatology.  This  method  is 
primarily  used  when  recent  observations  in  an  area  are  unavailable.  Next  up  the  ladder 
of  complexity  is  persistence.  This  method  is  used  when  recent  temperature  profiles  are 
available.  Uppermost  in  complexity  are  the  thermo-dynamic  models.  An  example  of  such 
a  model  is  the  Mellor-Yamada  Type  II  Model  which  is  used  in  FLENUMOCEANCEN’s 
(Fleet  Numerical  Ocean  Prediction  System)  operational  northern  hemisphere  ocean  fore¬ 
cast  model  TOPS  (Thermal  Ocean  Prediction  System). 

v 

During  this  past  performance  period^  we  ^examined  the  ability  of  statistical  techniques 
to  forecast  upper  ocean  deviations  from  climatology  on  1  to  7  day  time  scales,  deviations 

in  the  2.5  to  62.5  meter  region  in  depth  were  predicted  using  available  information  on 

•  '  v 

surface  conditions  and  “recent”  deviations  at  the  same  location.  The  results  represented 
the  best  predictions  that  regression  procedures  could  give  when  blending  surface  data  with 
recent  profiles.  Methods  were  developed  and  tested  on  data  collected  during  1983  and  1984 
from  OWS  Romeo  in  the  North  Atlantic. 

Information  on  changes  in  sea  surface  temperature  (SST)  could  be  combined  with 
recent  deviations  from  climatology  to  improve  on  persistence  in  the  upper  2.5  to  17.5 
meter  depths.  Below  17.5  m,  changes  in  SST  and  temperature  were  too  slightly  correlated 
for  linear  prediction  of  the  latter  to  be  effective.  Surface  forcing  information,  including 
wind  stress  and  heat  fluxes,  were  examined  for  their  usefulness  in  prediction.  While  wind 
stress  was  correlated  with  temperature  changes  near  the  surface,  the  ability  of  statistical 


predictions  which  incorporated  knowledge  of  wind  stress  was  little  better  than  persistence. 
To  improve  on  persistence,  we  believe  that  statistical  methods  require  information  as 
strongly  correlated  with  temperatures  below  the  surface  as  SST. 

We  also  used  a  one-dimensional  mixed-layer  model  which  employed  the  Mellor-Yamada 
Level  II  turbulence  closure  scheme  as  used  in  TOPS  to  simulate  shear  observed  in  the  upper 
ocean  during  MILE  (Mixed  Layer  Experiment).  In  the  upper  20  m  the  magnitude  of  the 
predicted  shear  agreed  fairly  well  with  that  observed.  There  was  also  some  agreement 
in  the  variability  of  the  shear  due  to  surface  forcing.  Between  20  and  30  m,  the  model- 
predicted  shear  increased  only  during  the  strong  wind  events  which  deepened  the  mixed 
layer  to  those  depths  while  the  observed  shear  had  a  more  constant  average  magnitude  and 
much  more  variability  on  short  (sub-inertial)  timescales.  Below  30  m,  which  was  below 
the  deepest  penetration  of  the  mixed  layer,  the  model  predicted  very  little  shear.  The 
observations,  however,  showed  considerable  shear  between  30  and  50  m  which  was  the 
depth  of  the  seasonal  thermocline.  The  average  magnitude  of  the  shear  in  this  region  was 
comparable  to  or  larger  than  that  in  the  mixed  layer  and  most  of  the  variability  was  on 
inertial  and  shorter  timescales.  A  spectral  analysis  of  shear  in  the  seasonal  thermocline 
showed  the  largest  peak  at  the  local  inertial  frequency  and  significant  contributions  from 
the  semi-diurnal  internal  tide  and  higher  frequency  motions. 

It  was  speculated  that  most  of  the  descrepency  between  the  observed  and  predicted 
shear  below  20  m  was  due  to  the  dispersion  of  inertial  motions  and  to  internal  waves. 
Between  20  and  30  m  the  inertial  motions  were  likely  generated  by  a  mixture  of  direct 
local  wind-driving  and  dispersion  with  dispersion  helping  to  maintain  the  level  of  the 
shear  when  the  mixed-layer  was  shallow.  Below  30  m  the  inertial  motions  must  be  caused 
primarily  by  dispersion  since  vertical  mixing  did  not  extend  to  these  depths  during  MILE. 

Using  24-hourly  output  of  temperature,  salinity  and  velocity  from  TOPS,  we  displayed 
the  performance  of  a  statistical  study  of  the  mean  and  variance  distributions  of  several 
upper  oceanic  environmental  parameters:  temperature,  temperature  variance,  mixed  layer 
depth,  Brunt-Vdisald  frequency  ( N ),  scalar  shear  (5),  Richardson  number  ( N2/S 2)  and 
eddy  diffusion  coefficient.  The  computations  covered  a  90-day  period  from  February  1  to 


May  1,  1985,  at  eight  ocean  stations.  Seven  of  these  locations  corresponded  to  weathership 
locations:  Papa,  November  and  Victor  in  the  Pacific,  and  Romeo,  Lima,  Charlie  and  Mike 
in  the  Atlantic. 

The  analysis  of  the  results  included  frequency  of  occurrence  distributions  over  various 
time  intervals  and  at  various  depth  levels  as  well  as  the  overall  distributions.  We  also 
examined  the  vertical  profiles  of  temperature  and  its  variance,  Brunt-Vaisala  frequency 
and  eddy  diffusion  as  a  function  of  month,  ocean  station  and  time  of  day. 

Finally,  we  used  a  hydrodynamic/ thermodynamic  sea  ice  model  to  predict  sea  ice 
concentration,  thickness  and  drift.  The  operational  model,  the  Polar  Ice  Prediction  System 
(PIPS),  is  based  on  the  Hibler  (1979)  model  only  with  the  addition  of  thermodynamics 
similar  to  that  used  by  Semtner(1976).  It  is  driven  by  atmospheric  forcing  from  the  Navy 
Operational  Global  Atmospheric  Prediction  System  (NOGAPS)  model,  as  well  as  monthly 

/'  .  j 

mean  climatological  ocean  currents  and  oceanic  heat  fluxes.  '  *tf  •'  " 

Prior  to  its  development,  empirical  models  using  the  interaction  between  ocean  cur¬ 
rents,  sea  ice  and  winds  were  used  at  FLENUMOCEANCEN.  The  first  operational  model 
was  that  of  Skiles  (1968).  It  related  ice  drift  to  geostrophic  wind  and  mean  upper  ocean 
currents.  It  was  eventually  replaced  by  a  model  by  Thorndike  and  Colony  (1982).  This 
linear  free  drift  model  used  an  updated  representation  of  ice  drift  based  on  numerous  ob¬ 
servations  of  drift,  geostrophic  winds  and  mean  ocean  currents.  However,  while  both  of 
these  models  predicted  ice  drift,  they  could  not  predict  ice  thickness  and  concentration. 
On  the  other  hand,  the  PIPS  model  couples  the  ice  movement  and  growth  to  thickness 
and  concentration  and  therefore  predicts  the  characteristics  of  the  ice  thickness  and  con¬ 
centration  as  well  as  temporal  changes  in  these  fields.  Results  using  the  model  are  given 
and  compared  to  observed  data. 
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I.  INTRODUCTION 

Traditional  approaches  (Jensen  and  Kuperman,  1982)  for  investigating  acoustic  wave 
propagation  in  an  ocean-like  waveguide  have  included  the  following  linear  theoretical  tools: 
rays,  normal  modes  (NM),  and  the  parabolic  equation  (PE)  approximation  to  the  linear 
wave  equation.  Ray  theory  contains  artificial  singularities  at  focii  and  caustics  and  does 
not  include  diffraction.  Normal  modes  become  cumbersome  when  extended  to  include 
range  dependence.  The  extensively  used  PE  method,  a  far  field  approximation  to  the  wave 
equation,  is  more  amenable  to  a  range  dependent  environment  (Tappert,  1977).  When  the 
propagating  signal  is  nonlinear  and/or  nonsinusiodal  in  time,  however,  the  NM  and  PE 
become  less  attractive  as  solution  techniques. 

An  alternative  to  these  methods  is  the  nonlinear  progressive  wave  equation  (NPE), 
which  has  been  shown  (McDonald  and  Kuperman,  1985;  McDonald  and  Kuperman,  1987) 
to  be  nonlinear  time  domain  equivalent  of  the  PE.  Specifically  formulated  to  study  weak 
shocks  and  finite  amplitude  pulses,  the  NPE  is  cast  in  a  frame  which  moves  at  a  constant 
speed  approximately  that  of  the  wave. 

The  NPE  is  derived  by  combining  fluid  momentum  and  mass  continuity  equations  and 
substituting  an  assumed  adiabatic  equation  of  state.  Retaining  lowest  order  nonlinearity 
and  assuming  propagation  within  a  narrow  cone  of  angular  half  width  8,  one  can  derive  a 
second  order  nonlinear  wave  equation  for  a  dimensionless  density  perturbation  R: 
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—  =  X72c2(R  +  0R2)  +  O(R3,R282),  (1) 

where  (3  —  1  4-  dlogc/dlogp ,  c2  =  dp/ dp,  and  R  =  p' / pQ;  c  is  the  local  sound  speed,  p  the 
pressure  and  0  is  a  thermodynamic  variable  whose  value  in  the  ocean  is  approximately  3.5. 
The  total  density  at  a  given  point  is  p  =  p0  +  p'  where  p0  is  constant  and  p'  is  a  small  but 
finite  perturbation. 


Equation  ( 1 )  is  cast  in  the  rest  frame  of  the  fluid.  By  transforming  to  a  frame  moving 
in  the  x  direction  at  a  speed  c0,  an  approximate  first  order  equation  (the  NPE)  for  forward 
going  waves  can  be  derived  (McDonald  and  Kuperman,  1987).  Taking  c  =  c0  -f-  cj,  where 
c i  is  a  small  environmental  fluctuation,  the  NPE  is  stated  as 


=  -T^'R+lTR2)-Ci 


The  operator  D/Dt  denotes  time  differentiation  in  the  moving  frame.  The  lower  limit  Xf 
of  integration  is  taken  to  be  in  the  quiescent  medium  ahead  of  the  wave.  On  the  right 
hand  side  of  Eq.  (2),  the  cj  term  is  refraction,  the  quadratic  term  is  nonlinear  steepening 
and  the  integral  term  is  diffraction  (implicitly  including  geometric  spreading). 


Equation  (2)  is  integrated  numerically  using  Crank-Nicholson  for  the  diffraction  term, 
and  flux  corrected  transport  (FCT)  methods  (Boris  and  Book,  1973;  Zalesak,  1979;  McDon¬ 
ald,  Ambrosiano  and  Zalesak,  1985;  McDonald  and  Ambrosiano,  1984)  for  the  x  derivative. 
The  environmental  sound  speed  fluctuation  C\  is  user  defined  and  can  be  range  and  depth 
dependent.  R  in  Eq.  (2)  is  related  to  pressure  using  the  approximate  equation  of  state  for 
water, 

2 

p=^[(/?+  l)n-ll,  where  n  =  7.15  (3) 

n 

As  a  linear  model  (0  =  0  in  Equation  (2)  ),  the  NPE  has  many  advantages  over  the 
frequency  domain  models.  Essentially  derived  as  a  wave  equation,  the  NPE  may  be  used 
to  find  the  solution  of  a  problem  with  a  wide  band  of  frequencies.  For  a  model  like  the 
parabolic  equation  (PE),  on  the  other  hand,  only  a  single  frequency  may  be  solved  for 
at  a  time  with  the  solutions  for  the  frequencies  superimposed.  The  PE  also  assumes  a 
continuous  sinusoidal  wave  as  the  source,  while  the  NPE  is  more  amenable  to  a  pulse. 
The  Normal  Modes  method  has  difficulty  with  a  range  dependent  environment,  and  Ray 
Theory  does  not  include  diffraction.  These  problems  do  not  exist  with  the  NPE.  Of  course, 
the  primary  objective  behind  the  NPE  is  to  solve  nonlinear  problems  for  a  long  and  narrow 
wave  guide. 


w 


£ 


y 


I 


/  ’ 


>• 

>;• 


V, 

V 

I* 


* 


£ 


r, 

K? 


II.  MODIFICATIONS  TO  CODE 


A.  Damping  Region 

To  prevent  the  reflection  of  sound  waves  from  the  bottom  of  the  simulation  grid,  a 
damping  region  is  needed.  Figures  1  and  2  show  two  possible  orientations  for  a  simulated 
ocean  environment,  the  former  for  a  shallow  ocean,  the  later  for  a  deep  ocean.  In  a 
shallow  ocean,  the  sediment  is  for  the  most  part  considered  infinite  and  therefore  does  not 
reflect  as  a  hard  surface.  In  a  deep  ocean,  the  ocean  is  deeper  than  the  simulation  grid 
extends,  eventually  reach  the  same  infinite  sediment  found  in  the  shallow  ocean  case.  In 
simulating  either  environment,  there  is  an  "artificial”  boundary  formed  by  the  bottom  of 
the  simulation  grid.  As  it  is  usually  undesirable  to  allow  reflections  from  this  boundary,  a 
damping  region  is  essential. 

To  implement  a  damping  region,  a  damping  term  must  be  added  to  the  right  hand 
side  of  Eq.  (2)  of  the  form 

(4) 

where  v(z)  is  a  depth  dependent  decay  factor 

v{z)  =  v0tanh(z  /  Xdz).  (5) 

Here,  vQ  is  the  decay  constant  and  Xdz  is  the  characteristic  decay  length.  To  prevent 
damping,  the  decay  constant  is  set  to  zero. 

B.  Range  Dependent  Sound  Speed  Profile 

In  shallow  ocean  simulations,  a  boundary  between  the  ocean  and  its  sediment  must 
be  defined.  One  way  of  including  a  sediment  into  a  simulation  is  to  use  a  range  depen¬ 
dent  sound  speed  profile  to  define  that  boundary  with  the  ocean  having  one  sound  speed, 
the  sediment  having  another.  We  have  therefore  altered  our  NPE  code  to  incorporate  a 
function  defining  the  depth  of  the  sediment  as  a  function  of  range.  This  function  is  called 
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each  timestep,  and  the  sound  speed  profile  of  the  simulation  region  is  adjusted  accordingly. 
With  the  simulation  frame  moving  at  a  sound  speed  c0,  the  array  cx  is  the  variation  of 
the  sound  speed  from  c0  at  each  point  on  the  grid.  As  the  frame  moves  forward  in  range, 
it  falls  upon  a  different  portion  of  the  ocean  and  sediment  being  simulated,  and  therefore 
the  sound  speed  profile  C\  is  changed. 

Once  the  array  cx  is  determined,  it  is  used  to  solve  for  refraction  in  the  NPE.  As  the 
sound  speed  variations  do  not  affect  diffraction,  we  may  exclude  the  integral  on  the  right 
hand  side  of  Eq.  (2).  The  equation  may  then  be  rewritten  as 


^  =  where  u  =  (cx  +  ~f?).  (6) 

This  equation  is  solved  using  a  one-dimensional  Flux  Corrected  Transport  scheme  (Boris 
and  Book,  1973;  McDonald  and  Ambrosiano,  1984;  Zalesak,  1979)  which  solves  for  in 
range,  i,  for  each  j  in  depth. 

C.  Diffusion  Timestep  Constraint 


The  NPE  uses  a  variable  timestep  with  restrictions  placed  on  it  by  both  the  Crank- 
Nicholson  integration  and  the  Flux  Corrected  Transport.  While  we  have  always  used  a 
timestep  constraint  for  the  nonlinear  FCT,  we  have  only  recently  added  one  for  diffraction. 
It  is  determined  by  excluding  the  two  left-most  terms  on  the  right  hand  side  of  Eq.  (2), 
leaving 


dR 

dt 


d 2 

dx—R. 

dz2 


(7) 


for  an  (x,z)  Cartesian  coordinate  system.  Writing  the  above  in  the  form  of  a  continuity 
equation  gives 


dR 

dt 


dz  l  2 


R 


which  has  the  general  form 


(8) 


%  =  -Tz(pv)'  where  v 


From  Eqs.  (8)  and  (9),  one  has 
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and  the  timestep  is  limited  by 


St  =  min(6t0 ,  St, 


CFL) 


where  St0  is  the  initial  timestep,  and  CFL  is  the  maximum  Courant-Friedrichs-Lewy  con¬ 
dition  set  by  the  user.  Though  the  Crank- Nicholson  integration  is  stable,  the  inclusion  of 
a  timestep  constraint  insures  greater  accuracy  as  St  is  adjusted  each  timestep  to  be  less 
than  a  predetermined  limit. 

For  completeness,  we  see  that  the  size  of  St  is  set  by  Flux  Correction  in  a  similar 
manner.  Excluding  the  diffraction  term  in  Eq.  (2)  gives 


dt  dxx  1  2  ’ 


Again  in  the  form  of  Eq.  (7),  one  may  write 


^mai  —  max ij(c\  {j  -(-  2  F-ij) 


for  all  points  i  and  j  on  the  grid  at  time  t. 


St  =  min(St0,St, 


D.  Tensile  Cutoff 


CFL). 


In  the  real  ocean,  the  pressure  at  any  given  point  is  never  less  than  zero.  In  a  sim¬ 
ulation,  however,  this  is  not  the  case.  Unless  restricted  not  to  do  so,  the  pressure  may 
become  negative.  To  prevent  this  from  happening,  a  relation  between  the  overdensity  and 


the  pressure  is  used  to  determine  when  the  overdensity  drops  to  a  level  below  which  the 
pressure  is  zero.  The  overdensity  is  then  set  to  the  overdensity  corresponding  to  a  zero 
pressure. 

Assuming  a  zero  pressure  for  an  undisturbed  fluid,  an  approximate  adiabatic  equation 
of  state  for  water  (the  equation  of  Tait)  is 

p,=  £2d[(A)"_1]  (is) 

Tl  p0 

where  n  —  7.15,  pi  is  the  pressure  of  the  disturbance,  p0  is  the  density  of  water,  p  is  the 
overdensity  of  the  disturbance  and  cQ  is  the  speed  of  sound  in  the  water  (Cole,  1948). 
Therefore,  the  total  pressure  is  p  =  pt  +  p0  where  p0  =  p0g(za  —  z)  is  the  undisturbed 
pressure,  g  is  the  acceleration  due  to  gravity,  za  is  the  depth  assigned  to  the  ocean  surface 
and  z  is  the  depth  of  the  given  point  at  which  the  pressure  is  being  determined.  Therefore 
a  zero  pressure  corresponds  to  an  overdensity  of 

Pmm  -  [l  -  ~z)]n  -  1  (16) 

co 

We  simply  set  the  overdensity  equal  to  pmin  when  the  overdensity  falls  below  pmtn. 

E.  Link  to  CSQ  Code  at  Sandia 

Sandia’s  CSQ  computer  code  is  used  to  solve  nonlinear  ocean  acoustics  problems,  yet 
unlike  the  NPE,  it  is  only  useful  for  short  range  simulations.  As  the  NPE  is  best  suited  for 
long  range  problems,  using  the  short  range  results  of  the  CSQ  is  a  useful  way  of  obtaining 
accurate  initial  density  perturbations  as  input  for  the  NPE.  The  grids  for  the  two  models 
may  not  coincide,  however.  We  therefore  use  a  bilinear  interpolation  in  both  range  and 
depth  to  interpolate  the  CSQ  grid  onto  the  NPE  grid.  Any  NPE  grid  points  falling  beyond 
an  edge  of  the  CSQ  grid  are  assigned  the  value  of  the  closest  edge  at  the  corresponding 
depth  (or  range). 

Output  from  the  CSQ  code  consists  of  the  arrays  RADC SQ(IC SQ),  ZCSQ(JCSQ) 
and  RCSQ(ICSQ,  JCSQ)  where  the  former  two  of  these  are  the  radial  distances  and 


depths  of  the  grid  points  respectively,  the  grid  being  ICSQ  points  in  range  by  JCSQ  points 
in  depth.  The  normalized  overdensities  at  the  grid  points  are  RCSQ(ICSQ,  JCSQ).  The 
required  input  for  the  NPE  code  (along  with  other  parameters  such  as  number  of  timesteps, 
size  of  initial  timestep,  etc.)  consists  of  RANGE,  NX,  N Z ,  DX  and  DZ.  These  are, 
respectively,  the  initial  range  of  the  simulation  grid,  the  number  of  x-grid  points,  the 
number  of  z-grid  points,  the  size  of  the  increment  in  range  and  the  size  of  the  increment  in 
depth.  Using  all  of  the  above  as  input,  we  interpolate  the  values  of  RCSQ(ICSQ,  JCSQ) 
onto  the  NPE  grid,  RN PE(NX,  NZ),  which  is  then  used  as  input  for  the  NPE. 
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III.  SIMULATIONS  OF  DIFFERENT  ENVIRONMENTS 


A.  Shallow  Ocean:  Phase  Shifting  at  Fluid-Fluid  Interface 

1.  Reflection  and  Transmission:  Linear  Theory 

Consider  a  fluid-fluid  interface  as  pictured  in  Figure  3  where  fluid  1  has  a  sound  speed 
of  co,  fluid  2  of  c'  =  c0  +  ci  and  Co  <  d .  Snell’s  law  gives  for  linear  plane  wave  propagation, 

Co  _  cos  a 
c'  costp 

Let  the  incident,  reflected,  and  transmitted  waves  be 


ei(k0(zcosa-zsina)-u>t)  (Jg) 

R(a)  ei(ko(zcoscl+z,ina)-u,t) ,  and 
T(tt)  gi(,k'(xco3ifi—zain<p)—wt) 

where  R(a)  and  T(o)  are  reflection  and  transmission  coefficients.  The  notation  conflict 
between  the  density  perturbation  R  and  the  reflection  coefficient  R(a)  is  resolved  by  the 
argument  a  for  the  reflection  coefficient.  Interface  conditions  are  continuity  of  pressure 
and  its  normal  derivative. 

Solving  for  the  reflection  coefficient  as  a  function  of  a,  we  find  that 


R(a)  = 


where 


'  <"0  2 
— r  —  cosia. 

ci 


Likewise  for  the  transmission  coefficient, 


For  the  critical  angle,  (18)  and  (20)  yield 
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cosac  =  — ,  <p  =  0  (22) 

c 

Equations  (19)  and  (20)  predict  positive  reflection  for  a  >  ac,  and  negative  reflection 
for  a  — ►  0.  i?(a)  is  complex  with  unit  modulus  for  ac  >  a  >  0,  with  R(ac)  =  1  and 
72(0)  =  —1.  The  reflected  wave  then  suffers  a  phase  change  ranging  from  0  at  a  =  ac 
to  180  degrees  for  a  — ►  0.  The  transmitted  wave  becomes  evanescent  in  z  for  a  <  ac,  so 
that  no  energy  radiates  deep  into  the  lower  medium.  (For  a  more  detailed  discussion  of 
reflection  and  transmission  at  a  fluid-fluid  interface,  see  Ingenito,  Ferris,  Kuperman  and 
Wolf,  1978.) 

2.  Results 
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The  results  presented  here  use  the  NPE  to  simulate  an  idealized  ocean  configuration 
from  a  frame  moving  at  the  water  sound  speed  Co  (Figure  4,  Region  1).  The  ocean  bot¬ 
tom  (Figure  4,  Region  2)  is  idealized  to  be  a  level,  homogeneous  fluid  with  sound  speed 
c0  +  c\ .  Beneath  the  bottom  a  damping  zone  (Region  3)  has  been  used  to  absorb  down¬ 
ward  propagating  waves  which  would  reflect  from  the  bottom  of  the  grid  and  reenter  the 
simulation. 

Figure  4  shows  density  perturbation  contours  used  for  the  initial  condition  R(x,  z,  0) 
of  the  simulation.  Azimuthal  symmetry  is  assumed.  Shown  is  a  spherical  wave  whose  radial 
profile  is  a  half  sine  wave  with  density  perturbation  R  (Eq.  (2))  ranging  from  0.0  to  0.1. 
A  density  perturbation  of  this  magnitude,  but  not  necessarily  of  this  shape,  could  result 
from  an  underwater  explosion.  Range  is  taken  in  the  x-direction  and  depth  is  taken  in  the 
z-direction.  The  grid  spacings  are  6x  =  2m  and  6z  =  4m.  Region  1  describes  the  ocean, 
Region  2  the  sedimentary  ocean  bottom  and  Region  3  the  damping  region  with  depths  of 
400  meters,  80  meters  and  120  meters  respectively.  The  width  of  the  simulation  box  is 
200  meters.  The  sound  speed  in  the  ocean  is  1500  m/s  and  that  in  both  the  sediment  and 
damping  region  is  1600  m/s.  The  initial  time  step  is  .01  seconds  and  is  thereafter  adjusted 
to  be  no  greater  than  half  of  that  allowed  by  the  Courant-Friedrichs-Lewy  condition. 
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From  these  initial  conditions,  both  linear  and  nonlinear  runs  were  made  (/3  —  0  and 
3.5  respectively).  For  the  linear  case,  Figure  5  shows  snapshots  at  times  when  the  grid  had 
traveled  225,  750,  1200  and  2400  meters.  For  the  nonlinear  case,  Figure  6  shows  them  at 
205,  730,  1200  and  2650  meters.  Although  the  contour  scalings  change  from  frame  to  frame, 
it  is  quite  clear  in  both  cases  that  as  the  range  increases,  the  reflection  from  the  fluid  -  fluid 
interface  goes  from  being  positive  to  negative  (solid  contours  are  positive  values,  dotted 
ones  are  negative).  While  in  the  linear  case  there  is  already  quite  a  significant  negative 
reflection  at  a  range  of  1200  meters,  this  is  not  so  in  the  nonlinear  case;  the  reflected  waves 
have  amplitudes  more  than  half  those  of  the  incident  waves  in  the  linear  case  while  they 
have  amplitudes  only  about  a  sixth  those  of  the  incident  waves  in  the  nonlinear  case. 

A  close  look  at  Figures  5  and  6  reveals  another  subtle  but  significant  difference:  the 
nonlinear  wave  travels  slightly  faster  than  the  linear  wave.  (Equation  (2)  reveals  a  nonlinear 
contribution  to  the  effective  local  sound  speed).  Therefore,  the  effective  sound  speed  jump 
from  Region  1  to  Region  2  is  not  quite  as  great  in  the  nonlinear  case  as  it  is  in  the  linear 
one.  This  increases  the  effective  transmission  coefficient  across  the  interface  for  small  angle 
incidence.  For  a  <  ac,e  in  Eq.  (20)  is  imaginary,  and  j|T||  increases  with  co,  if  c'  is  held 
fixed.  Figures  5  and  6  reveal  that  more  energy  is  in  fact  dumped  from  Region  1  to  Regions 
2  and  3  in  the  nonlinear  case  than  in  the  linear  case.  (Since  the  grazing  angle  is  subcritical, 
the  energy  does  not  radiate  deep  into  the  bottom,  but  resides  in  a  skin  layer  where  incident 
and  reflected  energy  fluxes  are  in  equilibrium).  It  seems  quite  possible  that  the  difference 
in  effective  sound  speed  profiles  is  the  explanation.  This  difference  also  means  that  the 
critical  angle  is  smaller  in  the  nonlinear  case,  and  therefore  the  waves  will  have  to  travel 
a  greater  distance  before  negative  reflections  emerge  (Plante,  Ambrosiano  and  McDonald, 
1987). 

B.  Deep  Ocean 

In  deep  ocean  simulations,  there  is  no  sediment  separating  the  ocean  and  the  damping 
region  (see  Figure  2).  Figure  7  shows  the  initial  waveform  used  in  our  simulation.  The 
initial  range  of  the  simulation  is  4194  m,  the  width  of  the  grid  1050  in  and  the  depth  of 


the  grid  6000  m.  The  lower  400  m  of  the  grid  is  a  damping  region.  (Refering  to  Figure  2, 
nx  =  700,  nz  =  750  and  noabsorb  =  700.)  The  grid  spacings  are  6x  =  1.5  m  and  6z  = 
8.0  m,  and  the  initial  timestep  is  .005  seconds.  The  speed  of  sound  in  the  ocean  is  1540 
m/s  at  the  surface,  1495  m/s  at  a  depth  of  1000  m  and  1560  m/s  at  depths  of  5500  m 
to  6000  m.  At  intermediate  depths,  the  sound  speeds  are  linearly  interpolated  from  these 
values.  The  source  is  a  single-cycle  sine  wave  of  wavelength  15  m  placed  at  a  depth  of  100 
m.  With  the  required  boundary  condition  being  that  of  a  pressure  release  surface  ( p  =  0), 
a  negative  image  source  is  also  placed  100  m  above  the  surface.  The  speed  of  the  moving 
grid  is  1510  m/s  with  respect  to  the  fixed  ocean  frame.  The  maximum  overdensity  (p1  / p0) 
of  the  initial  sine  wave  is  0.01.  For  the  linear  case,  ft  =  0.0,  and  for  the  nonlinear  case,  0 
=  3.5. 

Figure  8  shows  contour  plots  of  the  overdensities  on  the  grid  after  it  has  traveled 
19294,  41944,  64594  and  72144  m  for  the  linear  case,  and  Figure  9  shows  the  overdensities 
at  the  same  ranges  for  the  nonlinear  case.  It  is  immediately  apparent  that  the  pulse  in 
the  nonlinear  case  remains  less  confined  than  in  the  linear  case.  In  particular,  the  caustic 
which  has  formed  at  a  range  of  64594  m  is  quite  noticable  in  the  nonlinear  case  yet  is  barely 
seen  in  the  linear  case.  It  is  interesting  to  notice,  however,  that  in  the  linear  case,  at  the 
depth  of  the  source  and  in  the  center  of  the  grid  (about  five  tick  marks  from  either  side 
of  the  grid),  the  contours  are  fairly  concentrated  while  in  the  nonlinear  case,  the  contour 
lines  are  concentrated  both  before  and  after  rather  than  at  the  center  of  the  grid.  The 
pulse  as  a  whole  has  followed  the  same  basic  path  in  both  the  linear  and  nonlinear  cases, 
yet  the  local  density  perturbations  have  been  altered  considerably. 

C.  Ice  Surface 

Figures  10  and  11  show  contour  plots  for  a  simulation  in  which  an  ice  layer  was 
included  near  the  ocean  surface.  Similar  to  the  sedimentary  bottom  in  the  shallow  ocean 
simulation,  the  ice  layer  is  simply  a  portion  of  the  grid  in  which  the  sound  speed  is  greater 
than  that  in  the  ocean.  In  our  ice  simulation,  nx  —  100,  nz  =  150,  dx  =  dz  =  5  m  and  the 
initial  timestep  is  .005  seconds.  The  interface  between  the  ocean  and  the  ice  (delineated 


by  the  horizontal  lines  on  the  grids  in  Figures  10  and  11)  is  100  m  below  the  surface,  and 
the  damping  region  extends  from  a  depth  of  500  m  to  the  bottom  of  the  grid.  The  source 
is  a  single  cycle  sine  wave  with  a  half  angle  of  20  degrees  and  maximum  overdensity  of  0.9. 
Again,  /3  =  0  for  the  linear  case,  (3  =  3.5  for  the  nonlinear  case. 

As  in  the  shallow  ocean  simulations,  the  difference  in  sound  speed  between  the  two 
fluids  (in  this  case  ice  and  ocean)  is  greater  in  the  linear  case  than  in  the  nonlinear  case. 
Again,  this  is  due  to  the  additional  nonlinear  velocity  of  the  pulse  in  the  ocean  which  is 
eliminated  in  the  linear  case.  Therefore,  the  critical  angle  should  also  be  greater  in  the 
linear  simulation  than  in  the  nonlinear  simulation,  and  more  of  the  pulse  should  penetrate 
the  boundary  in  the  nonlinear  case.  Figures  10  and  11  readily  confirm  this  effect  as  more 
contours  can  be  found  in  the  ice  in  Figure  11  than  in  Figure  10.  Decreasing  the  sound 
speed  jump  between  the  ice  and  ocean  would  only  increase  the  penetration  of  sound  waves 
into  the  ice. 


CHECKS  AND  DIAGNOSTICS 
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A.  Dispersion  of  a  wave  packet 

1.  Normal  Modes 

Large  amplitude  waves,  such  as  those  resulting  from  an  underwater  explosion,  interact 
differently  with  the  ocean  bottom  than  small  amplitude  waves  because  of  the  nonlinear 
contribution  to  the  local  sound  speed.  Because  of  this  effect  we  also  expect  the  dispersion 
of  a  nonlinear  pulse  to  be  different  than  in  the  linear  case.  In  this  section  we  present 
numerical  results  using  the  NPE  model  designed  to  highlight  some  of  these  differences. 

A  simple  study  of  dispersion  can  be  performed  using  the  NPE  in  a  manner  reminiscent 
of  actual  ocean  experiments.  Again  we  consider  a  shallow  ocean  configuration  as  pictured 
in  Figure  1.  We  assume  in  this  case  an  ocean  depth  of  200  m  together  with  200  meters 
of  sediment,  the  lower  two-thirds  of  which  is  damped.  At  a  depth  of  50  m,  we  initialize 
a  three-cycle  sine-wave  packet  with  a  wavelength  consistent  with  a  frequency  of  50Hz  for 
the  assumed  sound  speed  in  the  water  column  and  limit  the  angular  spread  to  38  degrees. 
The  sound  speed  is  taken  to  be  1500  m/s  in  the  ocean,  1550  m/s  in  the  sediment  and 
damping  region.  The  initial  perturbation  is  then  allowed  to  propagate  downrange.  Data 
from  a  hypothetical  array  of  fifteen  hydrophones  placed  at  various  downrange  locations  is 
collected  during  the  run  as  a  collection  of  time  series.  We  ran  this  configuration  out  to  20 
km  for  the  linear  case  (0  =  0)  and  repeated  the  simulation  for  the  nonlinear  case  (0  =  3.5). 
An  initial  amplitude  R  =  p/po  =  0.2  was  chosen.  The  linear  calculation  is  completely 
insensitive  to  the  chosen  amplitude.  In  the  nonlinear  case  this  choice  represents  a  very 
powerful  perturbation  such  as  an  explosion. 

Linear  wave  theory  analysis  (Jensen  and  Kuperman,  1982)  predicts  that  at  all  ex¬ 
cept  very  short  ranges  propagation  occurs  in  discrete  normal  modes  of  the  shallow  water 
waveguide.  In  normal  mode  analysis,  solutions  <j>  to  the  wave  equation  for  a  harmonic 
point  source  of  a  single  frequency  (for  the  case  of  azimuthal  symmetry)  are  expanded  in  a 
set  of  complete  normal  modes 
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where  p(zo)  is  the  density  at  the  source  depth  zq  and  H^^kr)  is  the  Hankel  function  of 
the  first  kind.  The  normal  modes  un(z)  are  eigenfunctions  of  the  equation 


+  [k2(z)  -  k2n\un{z)  =  0  (24) 

In  Eq.  (24),  k(z)  is  the  wavenumber  and  k 2  is  the  eigenvalue.  The  waveguide  being 
dispersive,  these  normal  modes  propagate  at  different  speeds.  Thus  the  modes  will  separate 
as  the  pulse  propagates  downrange.  The  normal  modes  of  the  waveguide  may  be  calculated 
numerically  using  the  SACLANTCEN  normal  mode  propagation  model  SNAP  (Jensen 
and  Ferla,  1979).  Using  the  same  parameters  as  in  our  simulation,  normal  mode  analysis 
predicts  the  three  modes  depicted  in  Figure  12. 

Figure  13  shows  the  hydrophone  time-series  plotted  as  a  function  of  reduced  time  for 
simulation  data  taken  from  the  linear  case.  Reduced  time  is  defined  for  each  vertical  array 
at  a  given  range  as  the  the  actual  time  minus  the  time  at  which  the  pulse  first  arrives  at 
the  hydrophone  array.  Figure  13(a)  shows  the  readings  for  the  hydrophone  array  at  600  m 
downrange.  Further  downrange  at  20  km,  the  modes  have  separated  completely  as  can  be 
clearly  seen  in  Figure  13(b).  The  dashed  lines  in  Figure  13(b)  indicate  approximate  arrival 
times  of  the  normal  modes.  Figure  13(c)  shows  the  normalized  amplitudes  at  each  of  these 
arrival  times  as  a  function  of  hydrophone  depth.  We  can  see  that  the  general  appearance 
of  the  modes  is  consistent  with  the  predictions  of  normal  mode  analysis.  The  position  of 
the  maximum  amplitude  of  the  first  mode  as  well  as  positions  of  the  nodes  in  the  second 
and  third  modes  are  in  good  agreement  with  theory. 

The  hydrophone  data  in  the  nonlinear  case  reveal  distinct  differences  from  the  linear 
results.  Figure  14(a)  shows  these  results  for  a  range  of  600  m.  One  difference  that  is 
immediately  apparent  is  the  obvious  steepening  evident  in  the  nonlinear  case.  The  indi¬ 
vidual  waveforms  are  also  more  complex.  A  comparison  of  normal  modes  between  linear 
and  nonlinear  cases  is  not  completely  justifiable  since  decomposition  of  the  acoustic  signal 


into  normal  modes  is  strictly  valid  only  for  linear  waves.  However,  some  comparison  can 
be  made  between  the  two  cases  for  the  following  reasons.  First,  since  the  compressibility 
of  water  is  slight,  nonlinearities  are  necessarily  somewhat  weaker  than  they  would  be  in 
other  highly-nonlinear  problems.  Second,  as  we  saw  in  previous  section,  one  effect  of  non¬ 
linearity  on  the  fluid-fluid  interface  is  to  make  the  bottom  more  transparent.  Therefore 
the  initial  amplitude  of  the  perturbation  is  substantially  reduced  after  a  few  encounters 
with  the  water-sediment  boundary.  This  rapid  loss  of  wave  amplitude  to  the  sediment  is 
an  important  aspect  of  the  nonlinear  case  which  we  will  quantify  below.  As  the  wave  am¬ 
plitude  in  the  nonlinear  pulse  decays  its  iteraction  with  the  waveguide  becomes  essentially 
linear. 

With  these  points  in  mind  we  examine  the  hydrophone  data  from  the  simulation  of  the 
nonlinear  wave  shown  in  Figure  14(b)  at  20  km  downrange.  Although  the  results  evidently 
differ  from  Figure  13(b),  we  can  make  out  the  arrival  of  the  first  two  modes  as  indicated 
by  the  dashed  lines.  The  third  mode,  because  of  the  low  amplitude,  is  less  descemable  in 
the  complicated  wake  of  the  pulse.  Plotting  the  hydrophone  signal  amplitudes  of  first  two 
modes  along  the  dashed  lines  yields  Figure  14(c).  Again  there  is  relatively  good  agreement 
between  these  modes  and  the  eigemnodes  pictured  in  Figure  12.  An  exception  is  mode 
one  identified  in  the  first  three  dashed  lines  of  Figure  14(b)  and  plotted  in  14(c).  While 
the  mode  profile  at  the  third  location  agrees  very  well  with  the  linear  result,  the  first  two 
profiles  are  distinctly  different  exhibiting  an  amplitude  maximum  at  a  shallow  depth  of 
about  70  m  rather  than  the  deeper  110  m  predicted  by  normal  mode  analysis.  Since  these 
slices  in  depth  are  taken  near  the  front  of  the  propagating  signal,  this  departure  from  the 
linear  result  may  be  an  indication  of  the  change  in  the  eigenmode  caused  by  the  nonlinearly 
enhanced  local  sound  speed  near  the  interface  and  associated  losses  to  the  sediment. 

The  results  of  the  normal  mode  comparison  can  be  summarized  as  follows.  Agreement 
is  very  good  between  the  eigenmodes  predicted  by  linear  analysis  and  the  modes  that 
emerge  during  propagation  of  a  wave  packet  in  the  linear  NPE  calculation.  The  results 
show  that  the  NPE  model  faithfully  reproduces  linear  theoretic  results  in  the  linear  regime. 
In  the  nonlinear  regime,  the  concept  of  normal  modes  is  inappropriate  when  the  wave 


amplitude  is  large  and  so  nothing  can  be  said.  Time  scries  of  the  amplitude  sampled 
near  the  source  depth  at  a  relatively  short  distance  downrange  exhibit  the  steepening 
expected  of  a  weak  shock  wave.  As  the  nonlinear  pulse  propagates,  the  combined  effects 
of  bottom  transparency  and  the  usual  attenuation  with  range  reduces  the  amplitude.  The 
wave  thus  transitions  to  the  linear  regime,  becoming  less  steep.  Far  downrange,  dispersion 
separates  the  waveform  into  normal  modes  recognizable  from  linear  theory.  However,  the 
signature  of  the  linear  wave  is  distinct  from  that  of  the  nonlinear  case  which  contains  the 
effects  of  steepening  and  subsequent  unsteepening  as  well  as  loss  to  the  sediment.  The 
first  eigenmode  also  appears  to  be  distorted  near  the  front  of  the  pulse  probably  owing  to 
bottom  losses. 

2.  Frequency  Spectra 

In  the  previous  section  we  examined  normal  modes  of  the  waveguide  by  sampling  time 
series  from  vertical  arrays  of  points  at  several  locations  downrange.  We  may  also  sample 
time  series  along  a  long  horizontal  array  of  collection  points  at  the  source  depth.  Fourier 
analysis  of  these  time  series  yields  interesting  information  about  the  frequency  components 
of  the  wave. 

We  repeated  the  simulations  of  Section  1  and  Fourier  analyzed  the  time  series  into 
frequency  spectra.  Figure  15  shows  the  spectra  at  specific  locations  downrange  for  the 
linear  case.  The  salient  feature  of  these  snapshots,  ranging  from  600  m  to  20  km,  is  that 
the  spectrum  remains  peaked  near  the  initial  central  frequency  of  50Hz  for  the  duration  of 
the  simulation  with  less  being  found  at  higher  and  lower  frequencies  at  greater  ranges. 

Figure  16,  depicting  the  nonlinear  results,  presents  a  striking  contrast.  While  the 
predominant,  frequency  remains  near  50  Hz,  sidelobes  containing  other  frequencies  are 
plainly  evident.  At  relatively  close  range,  as  shown  in  Figure  16(a),  there  are  a  number 
of  higher  frequencies  seen  as  high  as  500  Hz  or  more.  This  is  the  stage  at  which  the  wave 
steepens  into  a  shock  giving  rise  to  high-frequency  components.  As  the  wave  continues  to 
propagate,  the  composition  of  frequencies  changes.  At  5  km  downrange  and  beyond,  low 
frequencies  emerge  and  contribute  substantially  to  the  spectrum. 


In  summary,  comparison  of  frequency  spectra  in  linear  and  nonlinear  cases  reveals 
plain  and  obvious  differences.  The  linear  wave  tends  to  lose  its  higher  and  lower  frequency 
components  as  it  propagates  downrange  with  the  vast  majority  of  the  spectrum  being 
focused  in  a  narrow  band  centered  at  50  Hz.  The  nonlinear  wave,  while  losing  its  higher 
frequency  components  initially  originating  from  the  steepness  of  the  initial  pulse,  continues 
to  shift  toward  lower  frequencies  as  it  moves  downrange,  though  it  also  remains  peaked  at 
50  Hz.  The  loss  of  the  higher  frequencies  indicates  the  unsteepening  of  the  wave. 

B.  Error  Plots 

To  determine  the  accuracy  of  the  NPE,  we  test  it  against  an  analytic  solution  calcu¬ 
lated  by  placing  multiple  images  both  above  the  ocean  surface  and  below  the  bottom.  The 
images  are  place  such  that  the  surface  is  a  pressure  release  (p  =  0  )  and  the  bottom  has  a 
zero  derivative  (~^  =  0)  at  all  times.  There  is  no  damping  region,  and  the  speed  of  sound 
is  1520  m/s  everywhere  so  that  there  is  no  sediment  or  ice.  The  grid  spacing  is  dx  =  0.7  m 
and  dz  —  1.75  m,  and  the  number  of  grid  points  are  nx  =  100  and  nz  —  400.  The  initial 
timestep  is  0.1  seconds.  We  ran  two  simulations,  the  first  with  a  gaussian  of  half  angle  72 
degrees  and  initial  range  of  140  m,  the  second  with  a  gaussian  of  half  angle  3  degrees  and 
initial  range  of  6893  m.  Both  sources  are  placed  280  m  below  the  surface. 

Figure  17  shows  horizontal  cross-sections  of  the  simulation  grid  at  ranges  from  0.1  m 
to  4.2  m.  The  dashed  lines  are  for  the  pressures  of  the  analytic  solution,  the  solid  lines 
are  for  the  pressures  calculated  with  the  NPE.  Notice  that  the  two  completely  overlap  in 
Figure  17(a)  because  the  analytic  solution  is  used  as  input  for  the  NPE.  Notice  in  Figure 
17(b),  however,  that  the  two  diverge  somewhat  only  to  re-converge  in  Figures  17(c)  and 
(d).  The  rms  error  for  this  run  is  shown  in  Figure  18,  the  error  calculated  for  the  entire 
grid  each  timestep  throughout  the  simulation  and  assigned  the  the  range  of  the  left  edge 
of  the  grid.  Notice  that  at  a  range  of  0.9  m,  the  error  is  approximately  7  %.  The  error 
peaks  at  a  range  of  about  500  m  and  decreases  to  a  level  of  less  than  3  %  at  a  range  of 
4100  m. 

Figure  19  shows  the  cross-sections  for  a  simulation  starting  at  a  range  of  6893  m, 
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corresponding  to  an  initial  half  angle  of  only  3  degrees,  and  extending  to  a  range  of  16718 
m.  Notice  that  in  Figures  19(a),  (b),  (c)  and  (d),  the  NPE  and  analytic  solutions  are  nearly 
equal.  Figure  20  shows  that  throughout  the  simulation,  the  error  is  always  less  than  1  %. 
It  is  obvious  that  the  NPE  is  much  more  accurate  for  small  angles  of  propagation  than 
large  angles,  as  expected  since  in  its  derivation,  approximations  are  made  which  decrease 
its  accuracy  for  high  angles. 


V.  CONCLUSIONS 


In  summary,  we  have  continued  to  modify  the  NPE  both  to  make  it  more  efficient  and 
to  make  it  more  general.  We  have  introduced  a  damping  region  to  prevent  reflections  from 
the  bottom  of  the  grid,  added  an  option  allowing  for  the  inclusion  of  a  range  dependent 
sound  speed  profile  for  shallow  ocean  problems,  modified  the  timestep  constraint  to  adjust 
not  only  for  refection  and  nonlinearity  but  also  for  diffusion  and  provided  an  ability  to 
link  the  NPE  to  the  shorter  range  nonlinear  CSQ  computer  code  from  Sandia.  We  have 
examined  both  linear  and  nonlinear  cases  of  shallow  and  deep  oceans  as  well  as  those 
of  an  ocean  with  a  layer  of  ice  at  its  surface  and  have  seen  that  there  are  appreciable 
differences  between  linear  and  nonlinear  simulations.  By  placing  arrays  of  hydrophones  in 
range  and  depth,  we  have  also  examined  the  structure  and  splitting  of  modes  as  the  pulse 
travels  downrange;  in  the  nonlinear  case,  we  have  seen  that  there  is  a  shift  in  the  spectrum 
from  high  to  low  freqencies  as  the  shock  unsteepens  while  in  the  linear  case,  the  spectrum 
remains  primarily  peaked  at  the  initial  dominant  frequency.  Our  error  plots  have  shown 
that  there  are  significant  reductions  in  error  as  the  angle  of  incidence  of  propagating  waves 
is  decreased. 
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TOPS  (Thermodynamic  Ocean  Prediction  System 

I.  INTRODUCTION 

It  is  an  important  goal  of  oceanographic  research  to  better  produce  accurate  tempera¬ 
ture  profiles  at  given  locations  in  the  ocean  system.  While  these  predictions  may  be  made 
at  various  levels  of  sophistication,  climatology  is  the  simplest  forecast  when  recent  obser¬ 
vations  in  the  area  are  unavailable.  When  these  observations  are  available,  persistence 
is  used.  The  most  sophisticated  and  data  intensive  methods  are  thermo-dynamic  models 
such  as  the  Mellor-Yamada  Type  II  Model  (see  Appendix  A)  embodied  in  the  Thermal 
Ocean  Prediction  System  (TOPS,  see  Clancy  and  Poliak,  1983),  the  current  Fleet  Numer¬ 
ical  Oceanography  Center  (FLENUMOCEANCEN)  northern  hemisphere  ocean  forecast 
model. 

We  begin  our  work  by  investigating  the  ability  of  linear  statistical  predictions,  a 
technique  intermediate  in  its  data  requirements  to  persistence  and  thermo-dynamic  models. 
We  attempt  to  produce  current  temperature  profiles  at  an  ocean  station  using  information 
about  current  surface  conditions  and  recent  profiles  at  the  same  station.  Profiles  are  for 
depths  of  2.5  to  62.5  meters  and  are  predicted  using  past  profiles  1  to  7  days  old.  We 
search  for  that  surface  information  which  enable  us  to  predict  profiles  below  the  surface 
with  less  error  than  persistence.  The  potential  of  statistical  predictions  are  determined  by 
our  ability  to  collect  this  surface  information  on  a  routine  basis  and  by  the  range  of  depths 
in  which  the  predictions  improve  on  persistence. 

Several  authors  report  success  with  statistical  techniques,  and  some  are  applying  re¬ 
gression  techniques  to  ocean-atmosphere  systems  (Davis,  1976  and  1978;  Barnett,  1984; 
Barnett  and  Hasselmann,  1979;  Hasselmann  and  Barnett,  1981).  In  all  of  this  work  both 
the  predicted  variables  and  the  predictands  are  averages  over  large  regions  and  times  of 
one  to  several  months.  Forecast  lags  are  taken  from  one  month  to  one  year. 

In  Section  II  on  Statistical  Predictions  of  Temperatures  in  the  Upper  Ocean,  we  focus 
our  attention  on  predictions  on  much  smaller  scales.  Our  prediction  is  for  profiles  at 


a  single  location  with  forecast  lags  of  days  rather  than  months.  Mesoscale  and  smaller 
events  become  a  more  important  part  of  the  variability  in  the  temperatures.  In  these 
scales,  persistence  becomes  the  method  against  which  all  others  must  be  judged. 

Of  course  it  is  also  important  to  compare  the  predictions  of  a  forecast  model  with 
actual  observations  whenever  possible  in  order  to  provide  validation  of  the  model  and  help 
determine  its  limitations.  A  number  of  investigations  show  the  ability  of  Ekman-type 
models  to  predict  the  response  of  the  surface  wind-driven  current  to  local  wind  forcing 
(Pollard  and  Millard,  1970;  Kundu,  1976;  Daddio  et.al.,  1978;  Wam-Vamas  et.al.,  1981; 
Martin,  1982).  However,  there  are  relatively  few  comparisons  of  model- predicted  upper- 
ocean  current  shear  with  observations.  This  is  partly  due  to  the  lack  of  data  available,  a 
result  of  the  difficulty  involved  in  measuring  current  profiles  near  the  surface  region  of  the 
ocean.  It  is  also  partly  due  to  the  lack  of  attention  received  in  simulating  or  forecasting 
shear  in  the  ocean. 

Several  mixed-layer  experiments  have  been  conducted  in  recent  years  for  which  sur¬ 
face  forcing  and  subsurface  measurements  of  temperature  and  velocity  are  available.  These 
include  the  Joint  Air-Sea  Interaction  Experiment  (JASIN),  the  Mixed  Layer  Experiment 
(MILE)  and  the  Storm  Transfer  and  Response  Experiment  (STREX).  In  Section  III  on 
Simulations  of  Different  Environments,  we  describe  a  comparison  of  model-predicted  cur¬ 
rents  and  shears  with  those  obtained  during  MILE. 

Finally,  in  Section  IV  on  Real-Time  Shear  Predictions,  we  present  part  of  a  larger 
effort  to  evaluate  the  capability  of  TOPS  to  predict  some  of  the  environmental  parameters 
known  to  :  (1)  affect  Ekman,  inertial  and  geostrophic  shear,  (2)  affect  the  propagation  of 
internal  waves  and  (3)  have  a  strong  correlation  with  the  level  of  internal  wave  activity. 
Interest  in  the  prediction  of  these  parameters  has  increased  with  the  recent  growth  in 
research  on  the  generation  and  propagation  of  internal  waves  in  the  upper  ocean. 

The  specific  purposes  of  the  current  study  are  twofold:  (1)  to  illustrate  the  real-time 
behavior  of  the  upper  ocean  at  eight  oceanic  stations  (four  in  the  Pacific,  four  in  the 
Atlantic,  Figure  31)  over  an  annual  period,  and  (2)  to  study  the  statistical  properties  of 


the  results  of  the  analysis  for  a  90-day  period 
of  February  1  -  May  1,  1985. 


II.  STATISTICAL  PREDICTIONS  OF  TEMPERATURES  IN  THE  UPPER 
OCEAN  ON  SYNOPTIC  TIME  SCALES 


A.  Data 

Ocean  Weather  Ship  Romeo  (47°jV,  17 °W)  provided  information  from  nightly  (22 
GMT)  XBT  casts.  The  data  set  contained  observations  from  May  4,  1983  to  August  1, 
1984.  The  temperature  profiles  were  interpolated  to  the  nine  depths  of  Fleet  Numerical 
Oceanography  Center’s  TOPS  grid  (2.5,  7.5,  12.5,  17.5,  25.0,  32.5,  40.0,  50.0  and  62.5 
meters).  Approximately  65  %  of  the  nights  had  casts  recorded. 

Observations  at  each  station  were  assigned  to  a  “spring”  (warming)  or  “summer” 
(warm)  season.  Spring  was  defined  as  beginning  4  weeks  prior  to  the  appearance  of  the 
seasonal  thermocline  and  ending  4  weeks  after  that  appearance.  Summer  was  defined  as 
beginning  with  the  appearance  of  the  seasonal  thermocline  and  ending  with  the  onset  of 
the  fall  cooling  (or  the  end  of  the  data  in  1984).  Two  summers  and  springs  were  observed, 
allowing  us  to  design  a  forecast  using  a  season  of  one  year  and  test  it  in  that  season  of  the 
other  year. 

Parsimonious  representations  of  the  deviations  from  climatology  were  obtained  by 
fitting  empirical  orthogonal  functions  (EOF’s)  to  the  spring  and  summer  deviations  from 
climatology.  Three  EOF’s  were  chosen  to  represent  the  deviations  for  Springs,  and  three 
were  chosen  for  Summers.  The  three  Spring  EOF’s  can  be  used  to  represent  95  %  of  the 
variability  about  0  in  the  deviations  from  climatology,  while  the  three  Summer  EOF’s  can 
represent  90  %.  Figure  21  shows  the  time  series  for  the  coefficients  of  the  EOF’s  during 
Spring  of  1983,  which  is  typical  of  the  four  seasonal  series.  Clearly,  coefficients  nearby  in 
time  are  highly  correlated,  so  highly  that  persistence  is  a  likely  means  of  forcasting  on  time 
scales  at  least  up  to  a  week. 

Table  1  displays  the  root  mean  square  error  as  a  function  of  depth  in  climatology  and 
one-day-old  persistence  for  the  two  spring  seasons.  Surprisingly,  the  precision  in  persistence 
does  not  increase  greatly  with  depth.  The  same  is  true  in  summer  months,  as  can  also 
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be  seen  in  Table  1.  Day-to-day  variations  are  the  same  or  smaller  at  the  very  top  as  at 
the  very  bottom,  though  62.5  meters  is  below  the  seasonal  thermocline  during  most  of  the 
summer  season.  There  are  much  larger  variations  in  the  summers  in  the  30  to  50  meter 
range  at  about  the  depth  of  the  seasonal  thermocline.  In  this  region,  large  differences  can 
occur  due  to  small  changes  in  the  shape  of  the  thermocline,  leading  to  large  variability. 

VVe  seek  surface  information  which  will  help  us  predict  the  errors  in  persistence.  Re¬ 
ports  on  surface  forcing,  including  sea  surface  temperature  (SST)  from  bucket  readings, 
wind  stress  (WS),  solar  radiation  and  latent  and  sensible  flux  were  made  eight  times  each 
day,  though  many  of  these  reports  are  missing.  This  data  was  condensed  by  averaging 
reports  nearby  in  time  as  summarized  in  Table  2.  Table  3  gives  the  correlations  observed 
by  season  between  the  daily  changes  in  temperature  at  four  depths  and  the  changes  ob¬ 
served  in  the  surface  variables.  Only  changes  in  SST  and  WS  show  large  correlations  with 
changes  in  temperatures,  and  those  correlations  are  limited  to  the  region  very  near  the 
surface  which  is  almost  always  in  the  mixed  layer. 

In  Appendix  B,  we  model  the  multivariate  series  of  the  EOF  coefficients  (which  repre¬ 
sent  deviations  from  climatology)  together  with  series  of  surface  data  as  a  simple  multivari¬ 
ate  time  series.  We  hope  to  use  this  model  to  incorporate  surface  information,  producing 
forecasts  superior  to  those  of  persistence. 

B.  Results 

From  the  correlations  in  Table  3,  it  is  apparent  that  SST  and  WS  have  the  best 
chance  of  yielding  good  statistical  forecasts.  Therefore  autoregressive  models  were  fit  as 
described  in  Section  A  for  three  choices  of  surface  information  ( 02(t )).  These  models  will 
be  indicated  by 

Stat  (SS):  O2H)  is  the  average  sea  suface  temperature  on 

day  t  as  described  in  Table  2. 

02(f)  is  average  wind  stress. 

02(f)  is  not  used.  No  surface  information  is  included. 


Stat  (WS): 
Stat  (NC): 


Once  the  model  is  fit,  forecasts  can  be  computed  given  any  combination  of  past  profiles 
and  past  or  current  surface  information.  For  illustration,  we  suppose  that  the  forecaster 
has  a  pair  of  XBT  profiles  at  times  t  —  1  and  t,  and  surface  information  at  times  t  and 
t  +  k.  The  goal  is  to  project  the  profile  below  the  surface  at  time  t  4-  k,  that  is,  to  forecast 
Oi(f4-fc)  —  ^i(0  given  Z\(t)  —  Oi(t)  —  0\{t  —  1)  and  02(t+k)  —  02(t).  This  will  be  referred 
to  as  a  ’’forecast”  at  lag  k,  though  it  uses  information  on  surface  conditions  at  time  t  +  k. 
The  methods  used  when  past  surface  information  is  used  (a  true  forecast)  are  identical. 
While  the  information  in  the  extra  XBT  at  time  i  —  1  is  not  necessary,  it  is  potentially 
useful  if  the  model  fitted  indicates  correlations  across  several  days  in  the  changes  in  EOF 


coefficients. 


Table  4  shows  the  RMSE  of  prediction  during  the  Spring  Seasons.  For  the  statisti¬ 
cal  forecasts,  the  errors  are  the  cross  validation  errors.  That  is,  forecasts  were  made  in 
Spring  1984  using  the  forecast  designed  in  Spring  1983  and  vise  versa.  Thus,  the  errors 
approximate  the  skill  of  the  forecasts,  or  their  ability  in  data  other  than  that  on  which  the 
method  was  designed.  For  a  given  lag,  all  forecasts  including  climatology  and  persistence 
are  computed  for  the  same  set  of  days  t  +  k  where  Zi(t)  was  available  and  02(t  4-  k)  and 
02(t)  were  available  for  02  both  WS  and  SST.  RMSE  in  climatology  fluctuates  slightly 
due  to  sampling  error.  For  simplification,  RMSE  is  reported  as  the  SQRT(average  squared 
error)  in  three  depth  groups:  S  (shallow)  2.5  to  17.5  m,  M  (mid)  25.0  to  40.0  m  and  D 
(deep)  50.0  to  62.5  m. 


Comparing  the  errors  in  the  statistical  forecasts  to  those  of  persistence,  we  see  that 
Stat(SS)  successfully  improves  on  persistence  in  the  shallow  depths,  particularly  at  lags 
greater  than  1.  Below  the  shallow  layers,  no  improvement  is  discemable.  Despite  the 
correlation  of  WS  with  changes  in  temperature  near  the  surface,  Stat(WS)  shows  no  more 
skill  than  persistence  or  Stat(NC). 


Table  5  repeats  the  analysis  of  Table  4  for  the  Summer  forecasts.  Summer  errors 
in  persistence  are  much  larger  than  those  for  Spring,  as  can  easily  be  seen  in  Table  2. 
Stat(SS)  shows  some  slight  improvement  over  persistence  in  the  shallow  layers,  especially 
at  longer  lags.  However,  forecasts  in  the  mid  and  deep  range  are  worse  than  persistence, 


indicating  that  the  correlations  observed  in  one  summer  between  SST  and  mid  to  deep 
temperatures  may  be  inapplicable  in  another  summer.  Neither  Stat(WS)  nor  Stat(NC) 
improve  on  persistence. 


The  large  errors  in  climatology  and  persistence  observed  in  the  summers  are  partially 
due  to  a  cold  core  eddy  moving  under  the  station  June  23,  1984  to  July  27,  1984.  This 
event  greatly  increased  the  variability  of  daily  changes  in  temperatures.  Table  6  gives  the 
RMSE  of  predictions  recalculated  omitting  the  eddy  period,  since  the  correlation  structure 
of  temperature  changes  is  different  from  that  of  non-eddy  periods.  While  RMSE  decreases 
compared  to  Table  5,  the  relationship  of  errors  from  persistence  and  errors  from  Stat(SS) 
are  nearly  the  same. 

The  action  of  the  statistical  forecasts  can  be  better  understood  by  examining  the  actual 
prediction  for  a  day.  Figure  22  shows  the  actual  profile  for  day  171  of  the  study  (June  20, 
1983).  Four  day  old  persistence,  in  the  form  of  the  profile  from  day  167,  is  also  shown. 
Between  167  and  171  there  was  an  increase  of  \.\°C  in  the  mixed  layer  temperature,  while 
the  temperatures  below  the  thermocline  remained  relatively  unchanged.  Incorporating 
knowledge  of  the  change  in  SST  in  the  statistical  forecast  yields  a  prediction  which  closely 
follows  the  true  profile  in  the  mixed  layer,  the  true  profile  in  the  seasonal  thermocline 
and  persistence  below  the  thermocline.  Typically,  the  statistical  forecasts  work  best  in 
this  kind  of  situation,  namely  one  with  a  large  change  in  mixed  layer  temperature  and 
relatively  small  change  below  the  seasonal  thermocline. 

Table  7  displays  the  MFPE  for  the  three  models  in  each  season.  An  investigator 
working  with  a  single  season  (as  in  Spring  1983)  would  consider  the  trio  of  numbers  for 
that  season  in  deciding  the  relative  value  of  different  surface  information.  In  every  case, 
the  investigator  would  correctly  decide  that  including  SST  in  the  models  would  yield  the 
best  forecasts  in  an  independent  data  set.  However,  in  three  of  four  seasons  he  would  also 
decide  that  including  WS  would  give  better  skill  than  using  no  surface  information.  In 
cross-validation,  we  see  that  the  improvement  due  to  using  WS  is  small  or  non-existent. 
This  result  is  due  to  large  differences  in  the  correlations  of  changes  in  WS  to  changes  in 
temperature  in  different  years.  Changes  in  SST,  on  the  other  hand,  maintain  about  the 


same  correlation  with  changes  below  the  surface. 

C.  Overview 

The  focus  of  this  effort  was  to  fit  a  statistical  model  of  the  daily  changes  in  temperature 
deviations  from  climatology  so  that  forecasts  could  be  computed  of  thermal  profiles  up  to 
seven  days  in  advance.  The  original  intention  was  to  model  the  process  of  changes  in 
deviations  without  surface  information,  namely  to  fit  Stat(NC).  Cross  validation  errors 
showed  this  model  to  be  no  better  than  persistence,  however.  The  major  reason  seems 
to  be  that  the  record  of  recent  thermal  profiles  does  not  contain  information  which  could 
help  one  predict  the  occasional  large  change  in  the  temperatures.  Since  these  events  are 
so  much  larger  in  magnitude  than  most  of  the  daily  changes,  these  rather  rare  occurrences 
will  contribute  most  to  the  correlations  fitted  in  any  one  season.  Since  these  large  changes 
could  either  be  positive  or  negative,  the  sign  of  the  observed  correlations  may  easily  be 
positive  in  one  season  and  negative  in  another.  As  an  illustration,  in  Summer  1983  the 
sample  correlation  was  negative  between  changes  in  the  first  EOF  coefficient  from  day  t  —  1 
to  t  and  day  t  to  t  +  1.  Hence  a  forecast  was  built  in  which  large  positive  changes  in  the 
coefficient  from  day  t  —  1  to  t  were  compensated  for  by  a  predicted  negative  change  from 
day  t  to  t  +  1.  The  sample  correlation  in  1984  was  positive  and  yielded  forecasts  that  a 
large  positive  change  from  t  —  1  to  t  would  be  followed  by  another  smaller  positive  change. 
These  types  of  inconsistencies  were  present  in  both  Springs  and  Summers  even  when  the 
eddy  period  was  removed  from  Summer  1984. 

The  inclusion  of  surface  information  should  provide  additional  data  with  which  to  pre¬ 
dict  the  daily  temperature  changes.  However,  the  correlations  actually  observed  between 
variables  such  as  SST,  WS  and  changes  were  surprisingly  small.  Partly,  this  must  be  due 
to  noise  in  the  measurement  of  the  surface  data.  For  instance,  reported  SST  is  actually 
the  average  of  two  bucket  readings  made  one  and  four  hours  before  the  XBT  reading. 
Therefore,  they  would  normally  differ  from  the  XBT’s  SST  by  as  much  as  ±0.4°C. 

More  difficult  to  model  are  the  effects  of  mixed  layer  depth  on  correlation.  When  25 
m  lie's  within  the  mixed  layer,  we  should  expect  changes  in  SST  or  WS  to  correlate  well 


with  changes  in  temperature  at  that  depth  and  design  forecasts  accordingly.  Similarly, 
when  the  mixed  layer  depth  is  very  shallow,  forecasts  should  not  use  surface  information 
in  the  forecast  below  the  thermocline.  When  data  from  more  years  becomes  available,  it 
sould  be  possible  to  fit  more  complex  models  which  explicitly  control  for  the  depth  of  the 
mixed  layer. 

This  work  indicates  that  there  is  not  enough  information  in  the  recent  past  of  the 
profiles  alone  to  build  predictions  which  improve  on  persistence.  However,  regression 
methods  can  blend  partial  information  on  changes  at  the  surface  with  information  from 
recent  profiles.  This  gives  greatly  improved  prediction  in  the  occasional  (but  important) 
case  of  large  sudden  changes  in  the  profile. 


III.  COMPARISON  OF  MODEL-PREDICTED  SHEAR  WITH 
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OBSERVATIONS  DURING  MILE 

A.  Model  Description 

The  model  used  here  is  a  one-dimensional,  differential  (primitive  equation)  mixed- 
layer  model  which  employs  conservation  equations  for  temperature,  salinity  and  the  two 
horizontal  components  of  momentum.  As  a  one-dimensional  model,  variability  is  provided 
for  in  time  and  in  the  vertical,  but  the  horizontal  derivatives  are  assumed  to  be  small 
and  are  neglected.  Vertical  mixing  is  parameterized  by  using  the  Mellor-Yamada  Level 
2  turbulence  closure  scheme  to  calculate  vertical  eddy  coefficients.  The  model  is  driven 
by  surface  fluxes  of  momentum,  heat  and  moisture.  This  model  is  used  in  TOPS,  and 
a  description  is  available  in  a  number  of  references  (Mellor  and  Durbin,  1975;  Clancy 
and  Martin,  1981;  Clancy  and  Pollack,  1983).  The  model  used  here  is  as  described  in 
these  references  except  that  the  turbulence  length  is  scaled  by  the  factor  0.2  rather  than 
0.1  (Equation  11,  Clancy  and  Pollack,  1983).  This  change  has  been  found  to  give  a 
slight  improvement  in  the  agreement  with  a  number  of  observation  sets,  and  it  will  be 
incorporated  into  TOPS  in  the  next  update  of  that  model. 

B.  MILE  Data 

MILE  was  conducted  near  Ocean  Station  PAPA  in  August  and  September,  1977, 
specifically  to  study  the  mixed  layer  and  provide  data  against  which  models  could  be 
tested.  Station  PAPA  is  located  at  50 °N,  145°IV  in  the  eastern  North  Pacific.  Data 
taken  during  MILE  included  3-hourly  meteorological  observations  which  could  be  used  to 
calculate  surface  forcing  for  an  ocean  model  and  subsurface  observations  of  temperature 
and  current  velocity  every  112.5  s.  The  subsurface  temperature  and  current  were  measured 
at  depths  from  5  to  50  m  at  3  m  intervals  (except  for  17  m)  and  at  70,  92,  125  and  175 
m.  The  depths  of  these  measurements  provide  sufficient  resolution  in  the  upper  50  m  to 
allow  the  calculation  of  vertical  temperature  gradients  and  shear.  Data  from  MILE  have 
been  analyzed  by  a  number  of  investigators  including  Davis  et  al.  (1981a, b),  Warn-Varnas 
et  al.  (1982),  Levine  et  al.  (1983)  and  Rubcnstein  and  Newman  (19S3). 
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C.  MILE  Simulation 


The  mixed-layer  model  was  initialized  from  the  observed  temperature  for  August.  19. 
the  beginning  of  the  MILE  study  (Figure  23).  The  temperature  profile  for  this  time  of  year 
was  notable  for  the  strong  seasonal  thermocline  that  existed  at  depths  between  30  and  50 
m.  Most  of  the  thermal  stratification  occurred  at  this  relatively  shallow  depth.  Velocity 
was  initialized  to  zero. 

Salinity  variability  was  neglected  in  this  study.  During  the  summer  in  this  region,  the 
stratification  of  the  upper  ocean  was  dominated  by  temperature.  The  model  salinity  was 
initialized  to  a  constant  33  parts  per  thousand  and  the  surface  moisture  fluxes  were  set  to 
zero. 


The  model  was  forced  by  surface  wind  stresses  and  heat  fluxes  calculated  from  the 
meteorological  data.  Penetration  of  solar  radiation  was  described  using  an  extinction 
profile  for  seawater  Optical  Type  II  (Jerlov,  1976).  Figure  24  shows  the  surface  wind 
stress  and  heat  flux  and  a  comparison  of  the  observed  and  model-predicted  sea-surface 
temperature  (SST)  and  mixed-layer  depth  (MLD).  Figure  25  shows  the  east  and  north 
components  of  the  horizontal  velocity  at  a  depth  of  S  m,  which  is  generally  within  the 
mixed  layer. 

The  SST  (actually  the  5  m  temperature,  since  this  is  the  depth  of  the  shallowest 
measurements)  shows  the  steady  warming  of  the  ocean  that  occurs  late  in  the  summer.  On 
days  when  the  winds  are  light,  the  steady  warming  is  modulated  by  a  diurnal  fluctuation  of 
about  0.2  C  which  is  caused  by  daytime  solar  heating  and  nighttime  cooling.  Two  storms 
occurred  during  MILE.  The  first  on  August  22-23,  caused  a  deepening  of  the  mixed  layer 
to  about  30  m  and  a  drop  in  SST  of  1.2  C.  The  second  on  Aug  31  -  September  1,  which 
was  slightly  weaker,  caused  mixing  to  about  26  m  and  a  SST  drop  of  0.4  C. 

The  model  reproduced  the  variability  of  the  mixed-layer  temperature  and  depth  fairly 
well.  The  major  discrepancy  was  that  the  model  underestimated  the  amount  of  cooling  tha  t 
occurred  during  the  August  22-23  storm.  Analysis  of  the  oceanic  heat  budget  indicated 
that  some  of  the  discrepancy  was  due  to  advective  cooling,  which  was  not  included  in  the 
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model  simulation.  As  a  result  of  underpredicting  the  cooling  on  August  22-23,  the  SST 
predicted  by  the  model  remained  biased  by  about  0.4  C  for  the  rest  of  the  experiment. 
For  a  simulation  initialized  on  August  25,  after  the  storm,  the  model  tracked  the  observed 
SST  very  well  (Figure  26). 

Inertial  motions  dominated  the  variability  of  both  the  observed  and  predicted  surface 
current  during  MILE  (Figure  25).  These  circular  motions,  which  were  an  effect  of  the 
earth’s  rotation,  were  generated  by  changes  in  the  magnitude  and  direction  of  the  wind. 
Wind  changes  that  were  in  phase  with  the  inertial  motions  in  the  mixed  layer  tended  to 
build  them  up,  and  changes  that  were  out  of  phase  caused  damping.  Horizontal  gradients 
in  the  inertial  motions  due  to  the  ocean  density  structure  and  the  variability  in  the  forc¬ 
ing  caused  them  to  disperse  both  horizontally  and  vertically.  Since  the  one-dimensional 
model  used  here  assumed  horizontal  uniformity,  these  processes  were  not  included.  The 
dissipation  of  inertial  motions  in  the  mixed  layer  by  these  processes  was  parameterized 
by  damping  the  inertial  motions  numerically  (using  an  implicit  treatment  of  the  Coriolis 
term)  over  about  a  10-day  timescale.  The  inertial  motions  observed  during  MILE  had  a 
typical  amplitude  of  20-30  cm/s  and  a  period,  which  depended  primarily  on  the  latitude, 
of  about  15.6  hrs.  Both  the  amplitude  and  phase  of  the  predicted  inertial  motions  in  the 
mixed  layer  agreed  fairly  well  with  the  observations. 

D.  Comparison  of  Observed  and  Predicted  Shear 

Figure  27  shows  observed  and  predicted  shear  profiles  at  specific  times,  and  Figure 
28  shows  time  series  of  shear  between  depths  that  correspond  to  the  locations  of  the 
current  meter  observations.  The  shear  between  11  and  20  m  is  shown  because  the  current 
record  at  14  m  is  very  short  (less  than  4  days).  Estimates  of  the  uppermost  reliable 
frequency  response  of  the  current  meters  are  2  to  8  eph  (Halpem  et  ah,  1981).  In  order 
to  eliminate  the  considerable  high  frequency  variability  and  to  facilitate  comparison  with 
the  model  predicted  shear,  the  time  series  arc  filtered  to  remove  fluctuations  with  period 
less  than  30  min.  The  high  frequency  variability  is  due  to  noise,  such  as  mooring  motions 
and  measurement  errors,  and  to  surface  and  internal  waves.  Figure  29  shows  a  sample 


comparison  of  filtered  and  unfiltered  time  series  of  shear  between  26  and  29  m. 


Values  of  shear  observed  in  the  upper  ocean  during  MILE  range  from  0.0  to  0.10  1/s, 
with  average  magnitudes  between  0.01  and  0.04  1/s  (based  on  a  current-meter  spacing  of 
3  m).  The  shear  between  5  and  11  m  (Figure  2Sa,b)  was  in  or  near  the  base  of  the  mixed 
layer  during  most  of  the  experiment.  The  magnitude  of  the  observed  and  predicted  shear 
at  these  depths  and  between  11  and  11  and  20  m  (which  were  generally  across  the  base  of 
the  mixed  layer  as  shown  in  Figure  28c)  agreed  fairly  well.  Some  agreement  could  be  seen 
in  the  variability  on  synoptic  and  shorter  timescales.  Since  the  shear  in  and  at  the  base  of 
the  mixed  layer  was  due  largely  to  directly  wind-driven  Ekman  and  inertial  currents,  the 
shear  in  this  region  tended  to  be  correlated  with  the  surface  forcing. 

Below  20  m  the  observed  and  predicted  shear  showed  increasingly  less  agreement. 
Between  20  and  30  m  the  predicted  shear  (Figure  28d-g)  showed  two  periods  with  strong 
peaks  coincident  with  the  two  storms  that  occurred  during  MILE  which  deepended  the 
mixed  layer  to  these  depths.  Following  their  initial  excitation,  the  predicted  currents 
dissipated  with  an  e-folding  timescale  of  about  10  days  which  was  imposed  by  the  model. 
The  observed  shear  between  20  and  30  m  showed  much  more  variability  than  the  predicted 
shear  on  short  (less  than  24  hr)  timescales.  This  variability  did  not  generally  appear  to 
be  well-correlated  with  the  surface  forcing,  although  the  August  31  -  September  1  storm 
showed  up  strongly  in  the  observed  shear. 

Below  30  m  the  model  showed  very  little  shear  since  the  mixed  layer  did  not  deepen 
below  30  m  and  the  model  had  no  means  to  transfer  momentum  to  these  depths  other 
than  by  vertical  mixing.  However,  the  observed  shear  profiles  (Figure  27)  showed  large 
peaks  in  the  seasonal  thermocline  between  30  and  50  m  depth.  The  time  series  of  shear 
(Figures  28h-l)  showed  the  average  magnitude  of  the  shear  in  the  seasonal  thermocline  to 
be  comparable  to  or  larger  than  that  in  the  mixed  layer.  The  magnitude  of  the  shear  at 
these  depths  was  relatively  steady  on  a  timescale  of  several  days,  but  showed  considerable 
variability  on  inertial  and  sub-inertial  timescales. 

The  discrepancy  between  the  predicted  and  observed  shear  in  the  upper  30  m  was 


due  to  a  number  of  factors  including  errors  in  modelling,  errors  in  the  forcing  and  neglect 
of  components  of  the  current  such  as  the  geostrophic  (pressure-driven)  flow  and  tidal 
motions.  In  addition,  it  was  speculated  that  some  of  the  discrepancy  between  the  observed 
and  predicted  shear  in  the  upper  30  m,  and  most  of  the  error  below  30  m,  was  due  to 
inertial  wave  dispersion  and  internal  waves. 

E.  Inertial  Wave  Dispersion 

Rubenstein  and  Newman  (1982)  investigated  the  shear  in  the  seasonal  thermocline 
during  MILE  and  found  the  largest  component  to  be  due  to  inertial  motions.  Since  these 
motions  were  below  the  mixed-layer,  they  were  not  driven  directly  by  the  local  wind.  It  was 
likely  that  the  inertial  energy  in  this  region  was  due  to  the  dispersion  of  inertial  motions 
from  the  surface  layer. 

The  dispersion  of  inertial  motions  was  caused  by  temporal  and  spacial  variation  of  the 
surface  forcing  and  by  horizontal  variation  in  the  ocean  thermal  structure.  As  a  result  of 
the  horizontal  pressure  gradients,  there  was  a  horizontal  divergence  of  the  surface  currents 
and  an  excitation  of  inertial  motions  below  the  mixed  layer.  Such  coupled  inertial-gravity 
wave  motions  are  referred  to  as  inertial-gravity  waves.  These  waves  propagated  inertial 
energy  both  horizontally  and  vertically. 

During  MILE,  the  shear  throughout  the  water  column  was  probably  affected  by  the 
dispersion  of  inertial  energy  to  some  degree.  In  the  upper  15-20  m  the  inertial  currents 
were  primarily  driven  by  the  local  winds  and  so  were  reproduced  fairly  well  by  the  model, 
although  some  discrepancies  were  apparent.  Between  20  and  30  m  the  inertial  motions  were 
likely  due  to  a  mixture  of  both  direct  local  wind-driving  and  dispersion,  with  dispersive 
effects  helping  to  maintain  shear  after  the  mixed-layer  had  shallowed.  Below  30  m  it  was 
speculated  that  the  inertial  motions  were  caused  primarily  by  dispersion  since  vertical 
mixing  did  not  extend  to  these  depths. 


F.  Internal  Waves 


There  was  significante  variability  at  frequencies  higher  than  the  local  inertial  frequency 
(15.6  hrs)  in  the  shear  times  series  between  30  and  50  m  depth  in  the  seasonal  thermocline 
(Figures  28  and  30).  Such  variability  also  occurred  at  shallower  depths  but  with  reduced 
amplitude.  Most  of  this  variability  was  likely  due  to  internal  waves  and  may  have  been 
caused  either  directly  by  the  propagation  of  internal  waves  or  indirectly  by  interaction  of 
the  internal  waves  with  the  inertial  currents. 

The  most  prominent  internal  waves  observed  during  MILE  were  those  of  semi-diurnal 
tidal  period  (12.4  hrs).  Internal  waves  with  this  period  were  common  in  the  ocean,  and  they 
were  referred  to  as  semi-diurnal  internal  tides.  It  was  speculated  that  they  were  generated 
by  the  interaction  of  the  barotropic  tide  with  topographic  features  such  as  seamounts  and 
shelf  slopes.  During  MILE,  the  semi-diurnal  internal  tide  had  a  typical  amplitude  (vertical 
isotherm  displacement)  of  10  to  15  m.  This  amplitude  maximum  occurred  at  the  depth 
where  the  statificaiton  was  strongest,  in  the  middle  of  the  seasonal  thermocline,  which  was 
the  location  of  the  node  of  the  lowest  order  internal  mode.  Spectra  of  shear  in  the  seasonal 
thermocline  during  MILE  exhibited  a  significant  amount  of  variability  at  the  frequency  of 
the  internal  tide  and  its  first  harmonic  (Figure  30). 

G.  Forecasting  Shear  Below  the  Mixed  Layer 

The  comparison  between  the  observed  and  predicted  shear  showed  that  there  was  a 
significant  amount  of  shear,  especially  below  the  mixed  layer,  that  could  not  be  accounted 
for  with  a  one-dimensional  mixed-layer  model.  This  raised  the  question  of  how  this  shear 
could  be  predicted. 

Inertial  motions  made  up  the  largest  component  of  the  shear  observed  below  the  mixed 
layer  during  MILE.  However,  explicit  forecasting  of  inertial  shear  caused  by  dispersion  of 
inertial  motions  was  difficult  for  several  reasons.  The  speed  at  which  inertial  energy  prop¬ 
agated  downward  from  the  surface  depended  on  small  scale  (less  than  100  km)  horizontal 
variability  of  the  surface  forcing  and  ocean  thermal  structure  which  was  usually  not  well 
known.  Furthermore,  because  the  vertical  propagation  speed  was  slow  (typical  values  were 
estimated  to  be  less  than  a  few  meters  per  day),  inertial  motions  may  have  been  affected  by 


advection,  mixing  and  dissipation  to  an  unknown  degree.  And  since  inertial-gravity  waves 
propagated  horizontally  as  well  as  vertically,  deep  inertial  motions  at  a  given  location  were 
due  in  part  to  non-local  effects.  Observations  have  been  shown  that  the  horizontal  coher¬ 
ence  and  the  coherence  between  surface  and  deep  inertial  motions  in  the  open  ocean  could 
fall  off  rapidly  with  depth  (Pollard,  1980).  Only  when  the  horizontal  scales  governing  the 
generation  and  dispersion  of  inertial  waves  were  very  well  known  (eg.  for  a  storm  with 
well-established  structure  such  as  a  front  or  hurricane  -  Geisler,  1970;  Price  1981,  1983) 
was  there  some  possibility  of  predicting  the  dispersion  of  inertial  motions  explicitly. 

Explicit  forecasting  of  high  frequency  internal  waves  will  likely  be  beyond  the  limits 
of  possibility  for  the  foreseeable  future.  Forecasting  low  frequency  internal  waves  such  as 
the  semi-diurnal  internal  seems  within  the  range  of  possibility  but  faces  problems  similar 
to  forecasting  the  dispersion  of  inertial  motions.  The  motions  are  not  generated  in  a 
sufficiently  direct  manner  to  be  easily  predicted. 

With  significant  practical  problems  surrounding  the  explicit  forecasting  of  shear  below 
the  mixed  layer,  the  high  correlation  between  the  magnitude  of  deep  shear  and  density 
stratification  suggested  the  use  of  statistical  methods  to  forecast  the  shear  in  the  region. 
Rubenstein  and  Newman  (1982)  showed  some  success  in  predicting  15-min  average  values 
of  shear  during  MILE  using  a  statistical  model  based  on  density  stratification  and  vertical 
spacing.  Other  parameters  could  be  investigated  for  use  in  a  general  statistical  model 
such  as  the  recent  history  of  local  winds,  horizontal  variability  of  density  and  forcing  or 
climatology  of  shear  in  the  region. 


IV.  REAL-TIME  SHEAR  PREDICTIONS 


A.  Comments 


Two  general  comments  must  precede  our  discussion  of  the  results.  First,  the  FLENU- 
MOCEANCEN  output  was  available  only  at  00:00Z  Greenwich  time,  so  that  the  Atlantic 
stations  always  fell  into  the  local  nighttime  (with  a  time  spread  from  0  to  2  a.m.)  and 
the  Pacific  stations  always  fell  into  the  local  daytime  (with  the  time  spread  from  11  a.m. 
to  2  p.m.  local  time).  This  made  the  comparison  of  the  Atlantic  stations  with  the  Pa¬ 
cific  stations  somewhat  more  difficult  in  that  the  diurnal  effects  needed  to  be  taken  into 
account. 


The  second  comment  concerns  the  vertical  resolution  of  the  TOPS  grid.  Several 
authors  have  shown  both  in  the  published  and  unpublished  technical  literature  that  the 
computation  of  BV  frequency  ( N ),  shear  (S')  and  Ri  =  N/S2  is  strongly  dependent  on 
the  vertical  interval  over  which  the  temperature  and  velocity  gradients  are  evaluated. 
Rubinstein  (1984)  shows  that  a  significant  correlation  exists  between  N2  and  S2.  This 
correlation  improves  with  increasing  Sz,  partly  resulting  from  the  relative  importance  of 
noise  decreasing  with  increasing  6z.  Considering  the  fall-off  of  both  temperature  and 
shear  spectra  with  increasing  wave-number,  we  effectively  filter  out  components  with  wave- 
number  greater  than  k  =  2i t/Sz  from  our  distribution  function  analysis. 


B.  Mean  Temperature  Profile  with  Depth 

Figures  31-33  display  the  monthly-average  temperature  profiles  at  the  four  North 
Pacific  stations  for  the  months  of  February,  March  and  April,  respectively.  Figures  34-36 
display  the  corresponding  results  for  the  four  North  Atlantic  stations.  Throughout  the  90- 
dav  period,  only  two  stations  show  the  existence  of  a  strong  seasonal  thermocline:  stations 
November  (30°iV,  140°JV)  and  Victor  (34°jV,  164 °E),  both  in  the  North  Pacific.  Weak 
seasonal  thermoclines  are  observed  at  stations  Papa  (50°iV,  145° IT)  in  the  Pacific  and 
Mike  (66° A\  2 ° E)  in  the  Atlantic.  The  others  all  show  a  more  or  less  homogeneous  layer 
flown  to  the  model  depth  of  500  in. 
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Only  at  Station  B  southwest  of  the  Bering  Sea  (50°JV,  165° W)  does  a  negative  tem¬ 
perature  gradient  occur,  and  it  is  accompanied  by  a  negative  halocline  which  maintains 
stability.  It  is  confined  to  an  ~  50  m  thick  layer  below  the  mixed  layer  in  February  but 
then  widens  to  ~  75  m  and  100  m  in  March  and  April,  respectively.  Below  that  lies  an 
almost  homogeneous  region  throughout  the  90-day  period. 

Seasonal  changes  that  affect  the  whole  profiles  down  to  500  m  occur  only  at  the 
Atlantic  stations  Charlie  (53 °iV,  36°JV)  and  Mike,  with  about  an  ~  .75 °C  warming  at 
Charlie  and  an  ~  .75 °C  cooling  at  Mike,  respectively. 

Stations  November  and  Papa  display  an  ~  1°C  cooling  from  February  to  March  and 
no  change  from  March  until  April.  Station  B  shows  an  ~  .75°C  cooling  from  February  to 
March  and  an  ~  .5°C  warming  from  March  to  April. 

C.  Mean  Temperature  Variance  Profiles 

Figures  37-39  display  the  monthly-averaged  vertical  profiles  of  the  temperature  vari¬ 
ance  in  the  four  North  Pacific  stations  for  the  months  of  February,  March  and  April, 
respectively.  Figures  40-42  display  the  corresponding  results  for  the  four  North  Atlantic 
stations. 

The  profiles  for  February  can  be  roughly  divided  into  four  groups:  (a)  the  variance 
is  dominated  by  a  large  peak  at  the  bottom  of  the  mixed  layer  (station  November)  or 
at  some  depth  below  the  mixed  layer  (station  Mike);  (b)  the  variance  profile  follows  the 
mean  temperature  profile,  being  constant  in  the  mixed  layer  and  continuously  dropping 
off  below  it  (stations  Papa,  Victor  and  Bering  Sea  in  the  Pacific,  and  station  Romeo  in  the 
Atlantic);  and  (c)  the  profiles  are  composed  of  slowly  increasing  or  decreasing  segments 
(stations  Charlie  and  Lima  in  the  Atlantic). 

The  situation  is  notably  different  in  March.  The  dominant  profile  is  now  one  where 
modest  peaks  near  the  bottom  of  the  mixed  layer  are  superimposed  upon  otherwise  con¬ 
stant  profiles.  Stations  Papa  and  November  in  the  Pacific  still  resemble  class  (b)  of  the 
February  profiles,  and  statior.  Lima  has  an  almost  constant  profile. 


In  April,  the  situation  again  changes  somewhat.  Now  station  Victor  exchanges  shape 
with  Papa  and  November  in  the  Pacific,  and  the  peak  near  the  bottom  of  the  mixed  layer 
increases  considerably  in  the  Bering  Sea.  In  the  Atlantic,  station  Romeo  goes  back  to  a 
class  (a)  and  station  Charlie  to  a  class  (b)  profile  shape.  At  Lima  the  variance  becomes 
totally  constant  with  depth,  and  at  Mike  it  increases  toward  the  surface. 

D.  Mixed  Layer  Depth  Variability 


We  must  first  digress  here  on  a  definition  of  mixed  layer  depth.  In  the  past,  definitions 
have  been  based  either  on  some  characteristic  of  the  profile  for  temperature,  density,  sound 
speed  or  Richardson  number.  The  most  common  criterion  adopted  in  the  numerous  mixed 
layer  simulation  studies  is  the  depth  at  which  the  temperature  falls  by  a  given  value  ST 
below  its  surface  value.  The  usual  value  of  ST  employed  by  NORDA  and  Navy  Postgrad¬ 
uate  School  (NPGS),  who  have  performed  the  most  seasonal  simulations  in  the  past,  is 
ST  =  .2 °C  (Garwood  and  Adamec,  1982;  Martin,  1984).  Occasionally,  values  of  ST  —  .1  °C 
are  used,  and  lately  FLENUMOCEANCEN  has  adopted  the  value  ST  =  .5 °C  in  their  eval¬ 
uation  of  TOPS  products  in  order  to  accommodate  certain  weaknesses  in  their  operational 
acoustic  models.  Other  criteria  have  utilized  the  depth  at  which  the  BV  frequency  reaches 
a  certian  value,  or  the  density  or  sound  velocity  falls  by  a  prescribed  value  below  its  respec¬ 
tive  surface  value  (Molinelli  et.al.,  1981).  And  finally,  the  depth  at  which  the  Richardson 
number  Ri  exceeds  .23  can  also  be  taken  as  the  MLD,  since  at  this  depth  little  shear  exists 
to  produce  turbulent  mixing.  In  this  paper  we  have  adopted  the  ST  =  .2 °C  criteria. 

Figures  43  and  44  illustrate  the  time  evolution  of  the  MLD  at  the  four  Pacific  and  the 
four  Atlantic  sites,  respectively,  according  to  the  ST  =  .2 °C  criterion.  Figures  45  and  46 
present  the  corresponding  results  according  to  the  Ri-number  criterion.  According  to  the 
temperature  criterion,  there  is  a  marked  “spring  transition”  at  almost  all  of  the  Pacific 
stations  where  the  MLD  suddenly  shallows  from  a  winter  to  a  summer  thickness  within 
the  interval  of  a  few  days  (Garwood  and  Adamec,  1982;  Martin,  1984).  The  transition 
appears  to  be  more  pronounced  (i.e.  to  take  place  over  a  shorter  period  of  time)  at  the  two 
stations  that  lie  on  the  50°  latitude  circle.  However,  in  the  Atlantic  this  shallowing  takes 
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place  over  a  much  longer  period.  At  stations  Charlie  and  Romeo,  with  latitudes  50±3°A^, 
respectively,  this  period  is  of  the  order  of  a  month,  from  about  March  15  to  about  April 
15;  at  station  Lima  (57 °N)  it  is  about  two  months,  from  about  March  15  to  May  15;  and 
at  station  Mike  (66° N)  it  is  an  almost  continuous  process  from  about  March  15  to  June  1. 

A  closer  examination  of  the  time  evolution  of  the  MLD,  as  defined  by  both  criteria, 
reveals  that,  in  fact,  there  are  several  transitions.  These  transitions  may  signify  separate 
events  or  stages  of  the  same  event  separated  by  various  time  intervals.  A  close  parallel 
examination  of  the  forcing  functions  is  needed  to  make  some  conclusions  as  to  the  exact 
causes  and  nature  of  these  transitions. 

Concentrating  on  the  events  that  result  in  sudden  shallowing  trends  of  the  MLD  that 
endure  for  weeks,  we  have  constructed  four  tables  giving  the  times  of  these  marked  layer 
depth  transitions  and  the  mean  layer  depth  between  transitions.  These  tables  have  so  far 
been  derived  from  visual,  semi-quantitative  observation  and  as  such  have  been  influenced 
by  subjective  judgments.  Hopefully,  such  procedures  can  be  automated  with  objective 
criteria  in  the  future. 

Tables  8  and  9  present  results  for  the  ST  =  .2 °C  criterion  and  Tables  10  and  11  for  the 
Ri-number  criterion.  We  have  grouped  the  transition  periods,  ti,  to  fall  within  two- week 
periods  (i.e.  ti  ±  7  days).  Many  of  the  groups  have  a  much  smaller  spread  about  the 
mean.  Up  to  seven  such  periods  have  been  detected  in  the  Atlantic  and  six  in  the  Pacific 
(according  to  the  Ri  criterion). 

Examining  the  results  in  Table  12  (with  the  ST  criterion)  for  the  Pacific,  we  note 
immediately  that  the  March  transition,  so  dominant  in  the  Atlantic  is  absent,  and  only 
one  station  (Victor)  shows  a  transition  in  April.  The  May  transition  is  the  dominant  one 
and  occurs  at  all  four  sites,  with  f4  =  May  16.5  (±  6  days).  There  are  two  June  transitions, 
with  June  4  (±  3  days)  occurring  at  three  stations  and  June  18  occurring  only  at  November. 
It  is  interesting  to  note  that  the  two  station  pairs  located  on  50 °N  and  32 °N  (±2)  have 
their  transition  periods  each  5  days  apart,  respectively.  These  stations  constitute  the 
vertices  of  a  spherical  rectangle  since  the  latitudes  also  form  pairs  as  164.5°H/  (±.5)  and 


142.5° \V  (±2.5),  respectively.  On  the  perimeter  of  this  rectangle,  the  spring  transition 
moves  around  in  a  counterclockwise  fashion  with  respective  arrival  times  of  5,  3  and  5 
days  after  leaving  Victor.  The  stations  are  too  far  apart  to  ascribe  the  effect  to  synoptic 
weather  causes,  but  dynamic  ocean  effects  in  the  form  of  wave  propagation  may  possibly 
be  included.  Most  likely,  this  phenomena  is  only  accidental  and  the  transition  times  are 
determined  by  the  local  evolution  of  the  atmospheric  heat  fluxes. 

E.  Frequency  Distribution  of  Various  Parameters 

Statistical  frequency  distributions  of  mixed  layer  depth  are  shown  in  Figures  18-21 
for  the  months  of  February  through  April.  Two  different  kinds  of  mixed  layer  depths  are 
shown.  One  is  the  temperature  difference  criteria  of  definition.  This  definition  reflects  the 
mixing  signature  on  the  temperature  profile.  The  other  definition,  the  critical  Richardson 
number,  reflects  the  existing  depths  of  turbulence  penetration  which  in  turn  reflects  storm 
intensity  and  convective  cooling.  Figure  47  shows  a  mixed  layer  distribution  of  around  100 
meters  in  the  Pacific  with  the  exception  of  station  Victor  which  shows  a  distribution  of 
around  150  meters.  The  difference  at  station  Victor  can  be  attributed  to  its  close  proximity 
to  the  Ivuroshio  where  passing  eddies  tend  to  homogenize  the  temperature  structure. 

Figure  48  shows  the  analogous  situation  in  the  Atlantic.  The  mixed  layer  depths  are 
a  much  deeper  150  to  400  meters,  reflecting  the  weaker  stratification  that  exists  in  the 
Atlantic.  Figure  49  shows  a  turbulence  penetration  of  100  to  150  meters  in  the  Pacific 
which  agrees  with  the  maximum  mixing  signature  of  the  temperature  profiles  in  Figure 
47.  Figure  50  shows  the  analogous  situation  in  the  Atlantic.  The  gaps  in  the  critical 
Richardson  mixed  layer  depth  histograms  relect  the  absence  of  model  grid  points.  The 
grid  is  stretched  and  varies  in  resolution  from  a  few  meters  to  a  hundred  meters  in  deeper 
regions.  The  maximum  penetration  of  turbulence  agrees  with  the  maximum  mixing  depth 
signature  on  the  temperature  profiles  (taking  into  account  this  coarse  grid  resolution  at 
depth). 

Figures  51-56  show  profiles  of  the  Brunt- Vdisdld  frequency  versus  depth.  As  expected 
the  Brunt- Vdisald  frequency  peaks  in  the  thermocline  region  where  the  sharpest  tempera- 
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ture  gradients  exist.  First,  consider  the  Pacific  stations  for  the  time  sequence  of  February, 
March  and  April.  As  evident  from  Figure  43,  at  the  Bering  Sea  and  Victor  stations,  the 
spring  transition  begins  in  April.  This  shows  up  at  station  Victor  in  Figures  51-53  as  an 
eroding  of  a  deep  mixed  layer  and  thermocline  structure  with  the  Brunt- Vdisald  frequency 
showing  a  peak  at  what  used  to  be  the  mixed  layer  region.  At  the  Bering  Sea  station 
there  is  a  less  pronounced  change  in  the  mixed  layer  region.  At  stations  Papa  and  Novem¬ 
ber  there  is  no  noticeable  change  since  the  spring  transition  occurs  in  May,  as  shown  in 
Figure  43.  The  Atlantic  Ocean  analogue  is  depicted  in  Figures  54-56.  The  peaks  of  the 
Brunt-Vdisala  frequency  are  not  as  sharp  in  the  thermocline  because  the  stratification  is 
weaker.  Figure  44  suggests  that  we  should  see  the  effects  of  a  spring  transition  at  three  of 
the  stations.  We  find  this  to  be  the  case. 

Figures  57-64  show  shear  histograms  (frequency  distributions)  at  various  depths  and 
integrated  configurations  over  depth.  The  Pacific  histograms  show  higher  values  of  shear 
than  the  Atlantic  histograms.  The  histograms  also  show  that  higher  values  of  shear  tend 
to  occur  near  the  surface. 

Figures  65-70  show  histograms  of  the  Richardson  number  distributions  for  the  Pacific 
and  Atlantic  oceans.  In  the  Pacific  we  note  a  more  symmetric  distribution  than  in  the 
Atlantic,  reflecting  the  weaker  stratification  of  the  Atlantic. 

Figures  71-74  show  histograms  of  eddy  coefficients.  We  observe  that  the  coefficients 
are  larger  in  the  Atlantic  than  in  the  Pacific.  This  is  because  the  weaker  stratification  of 
the  Atlantic  enables  a  larger  turbulent  intensity  to  build  up,  which  in  turn  is  reflected  in 
a  larger  eddy  coefficient. 
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V.  CONCLUSIONS 


The  focus  of  our  study  in  Section  II  was  to  fit  a  statistical  model  of  the  daily  changes  in 
temperature  deviations  from  climatology  such  that  forcasts  could  be  computed  of  thermal 
profiles  up  to  seven  days  in  advance.  Although  our  original  intent  was  to  model  only 
the  process  of  changes  in  deviations,  ignoring  surface  information,  we  found  that  cross 
validation  errors  proved  this  model  to  be  no  better  than  persistence.  The  major  reason 
for  this  appeared  to  be  that  the  record  of  recent  thermal  profiles  did  not  contain  sufficient 
information  to  predict  the  occasional  large  change  in  temperatures.  The  inclusion  of  surface 
information  should  provide  additional  data  with  which  to  predict  the  daily  temperature 
changes,  yet  as  of  now  there  appears  not  to  be  enough  information  in  the  recent  past  of 
the  profiles  alone  to  build  predictions  which  improve  on  persistence. 

In  Section  III,  a  one-dimensional  mixed-layer  model,  which  employed  the  Mellor- 
Yamada  Level  2  turbulence  close  scheme,  was  used  to  simulate  the  upper  ocean  during 
MILE.  The  model  was  initialized  from  the  observed  temperature  and  forced  using  sur¬ 
face  wind  stresses  and  heat  fluxes  calculated  from  local  meteorological  observations.  The 
variability  of  the  mixed-layer  temperature  and  depth  were  predicted  fairly  well,  the  ma¬ 
jor  discrepancy  being  an  under-prediction  of  the  cooling  during  the  first  storm  in  which 
advection  appeared  to  play  a  role.  The  amplitude  and  phase  of  the  surface  wind-driven 
current  also  was  in  good  agreement  with  the  observations. 

A  comparison  of  observed  and  predicted  shear  showed  that  in  the  upper  20  m  (which 
contained  the  mixed  layer  during  most  of  the  experiment)  the  magnitudes  agreed  well. 
Between  20  and  30  m,  the  predicted  shear  increased  only  during  the  strong  wind  events 
which  deepened  the  mixed  layer  to  those  depths,  whereas  the  observed  shear  had  a  more 
constant  average  magnitude  and  much  more  variability  on  inertial  and  shorter  timescales. 
Below  the  deepest  penetration  of  the  mixed  layer  (about  30  m)  the  model  predicted  very 
little  shear.  The  observations,  however,  showed  considerable  shear  between  30  and  50  m, 
the  location  of  the  seasonal  thermocline  and  region  of  strong  thermal  stratification.  The 
upper  20  ii,  or  so  were  reproduced  fairly  well  due  to  the  inertial  currents  being  driven 
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mainly  by  local  winds.  The  discrepency  between  the  observed  and  predicted  shear  below 
20  m  was  probably  due  primarily  to  the  horizontal  and  vertical  dispersion  of  inertial  energy 
and  to  internal  waves.  Between  20  and  30  in  the  inertial  motions  were  likely  generated 
by  a  mixture  of  direct  local  wind-driving  and  dispersion  with  dispersive  effects  helping  to 
maintain  shear  when  the  mixed  layer  was  shallow.  Below  30  m  the  inertial  motions  were 
probably  caused  by  dispersion  since  vertical  mixing  would  not  extend  to  those  depths. 

Finally,  in  Section  IV,  a  data  report  on  real-time  shear  prediction  was  presented.  The 
oceanic  forecast  model,  TOPS,  generated  its  own  Ekman  and  inertial  velocity  and  received 
as  input  a  geostrophic  velocity  from  climatology  or  analysis.  We  forecasted  the  mean  and 
variance  distribution  of  temperature,  mixed  layer  depth,  Brunt-Vaisala  frequency,  scalar 
shear,  Richardson  number  and  eddy  diffusion  coefficients  for  the  period  of  February  1  to 
May  1,  1985.  The  forecast  was  performed  at  four  locations  in  the  North  Atlantic  and 
four  in  the  North  Pacific.  The  forecasted  mean  and  variance  distribution  could  be  used  to 
describe  the  shear,  acoustic,  internal  wave  and  fine-structure  oceanic  environment. 


HIBLER  ICE  MODEL 


I.  INTRODUCTION 


A  hydrodynamic/ thermodynamic  sea  ice  model  is  presently  being  used  at  the  Fleet 
Numerical  Oceanography  Center  (FNOC)  to  predict  sea  ice  concentration,  thickness  and 
drift.  The  operational  model,  the  Polar  Ice  Prediction  System  (PIPS),  is  driven  by  atmo¬ 
spheric  forcing  from  the  Navy  Operational  Global  Atmospheric  Prediction  System  (NO¬ 
GAPS)  model,  as  well  as  monthly  mean  climatological  ocean  currents  and  oceanic  heat 


fluxes. 


Prior  to  the  development  of  the  PIPS  model,  operational  Arctic  sea  ice  prediction 
at  FNOC  consisted  of  empirical  models  that  made  use  of  the  interaction  between  ocean 


currents,  sea  ice  and  winds.  The  first  operational  model,  that  of  Skiles  (1968),  related  ice 
drift  to  geostrophic  wind  and  mean  upper  ocean  currents.  This  model  was  replaced  by 
the  Thorndike  and  Colony  (1982)  model,  a  linear  free  drift  model  that  used  an  updated 
representation  of  ice  drift  based  on  numerous  observations  of  drift,  geostrophic  winds  and 
mean  ocean  currents.  Both  of  these  models  predicted  only  ice  drift;  no  information  on 
ice  thickness  or  concentration  could  be  obtained  from  them.  The  PIPS  model,  however, 
couples  the  ice  movement  and  growth  to  thickness  and  concentration.  As  a  result,  the 
PIPS  model  is  able  to  predict  the  characteristics  of  the  ice  thickness  and  concentration  as 
well  as  temporal  changes  in  these  fields. 


II.  FORCASTING  ICE  THICKNESS  AND  CONCENTRATION 

IN  THE  ARCTIC 


A.  Model  Description 

The  PIPS  model  is  based  on  the  Hibler  (1979)  model  with  the  addition  of  thermo¬ 
dynamics  similar  to  that  used  by  Semtner  (1976).  Five  main  components  are  used  to 
define  the  sea  ice  model:  the  momentum  balance,  ice  rheology,  ice  thickness  distribution, 
ice  strength  and  an  atmosphere-ice-ocean  heat  budget.  The  momentum  balance  used  to 
describe  ice  drift  is 


m 


Du 

Dt 


=  m/k  x  u  +  ra  -f  tw  —  mgVH  +  F, 


(25) 


where  m  is  the  ice  mass  per  unit  area,  u  is  the  ice  velocity,  /  is  the  Coriolis  parameter, 
ra  and  tw  are  air  and  water  stresses,  g  is  the  acceleration  of  gravity,  H  is  the  sea  surface 
dynamic  height  and  F  is  the  force  due  to  variation  of  the  internal  ice  stress.  Ice  is  considered 
to  move  in  a  two  dimensional  field  with  forcing  applied  through  simple  planetary  boundary 
layers. 


The  air  and  water  stresses  are  defined  as 


ra=PaCa\\Jg\lJ9  (26) 

Tw  =  pwCw\Uw  -  u|[(Uw  -  u)cos0  +  k(Uw  -  u )sin6]  (27) 

where  u  is  the  ice  drift  velocity,  U?  is  the  surface  wind  from  the  NOGAPS  model,  Uw  is 
the  geostrophic  ocean  current,  Ca  and  Cw  are  the  air  and  water  drag  coefficients,  pa  and 
pw  are  the  water  densities  and  6  is  the  turning  angle  for  the  water  (McPhee,  1975). 

The  ice  rheology,  a  viscous-plastic  constitutive  law,  relates  the  ice  stress  to  ice  defor¬ 
mation  and  strength  via  the  following  equation: 

°ij  =  2 „(£,j,P)£i]  +  [C(eij,^)  ~  r](el],P)]£kkSij  -  P6tJ/ 2  (28) 
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In  Eq.  (28),  <J{j  is  the  two  dimensional  stress  tensor,  £{j  is  the  strain  tensor,  P/2  is  a 
pressure  or  strength  term,  and  ?;  arc  nonlinear  bulk  and  shear  viscosities  and  is  a 
delta  function.  Ice  flows  plastically  for  normal  strain  rates  and  deforms  in  a  linear  viscous 
manner  for  small  strain  rates. 

Ice  thickness  evolves  as  a  result  of  both  dynamic  and  thermodynamic  effects.  The 
PIPS  model  uses  a  two  level  ice  thickness  approach.  Ice  is  either  thick  or  thin  with  the 
division  between  the  two  being  0.5  m.  The  compactness,  A,  is  the  area  within  a  grid  cell 
covered  by  thick  ice  while  (1  —  A)  is  the  area  covered  by  thin  ice.  Thickness  and  compaction 
are  defined  by 


d(uh )  d(vh) 


+  Sh  +  diffusion 


dA  d(uA )  d(vA) 


+  Sa  +  di  f  fusion 


where  Sh  and  Sa  are  thermodynamic  terms  defined  by 


Sh=f(-)A  +  (l-A)m 


(m{l-A)  +  P,  if  /(0)  >  0; 
\/3/  if  /(0)  <  0. 


In  Eq.  (32),  /?  is  defined  as 


"ft 


if  Sh  >  0; 

■ )Sh ,  if  Sh  <  0. 


with  f{h)  the  growth  rate  of  ice  thickness  h  and  hQ  the  0.5  m  demarcation  between  thick 
and  thin  ice. 

Ice  strength,  a  function  of  the  ice  thickness  distribution,  is  given  by 


P  =  P*hexp[-C(  1  -  A)}, 


(34) 


where  P*  and  C  are  fixed  empirical  constants. 

The  PIPS  model  also  has  the  ability  to  determine  the  growth  and  decay  rates  of 
ice  thickness.  The  vertical  growth  rates  of  both  thick  and  thin  ice  are  determined  by 
calculating  a  heat  budget  at  the  top  and  bottom  surfaces  of  the  ice  and  by  adding  the 
heat  absorbed  by  leads  via  lateral  melting.  Similar  to  the  formulation  of  Semtner  (1976), 
heat  is  transferred  through  the  ice  by  assuming  a  linear  temperature  profile  along  with  a 
constant  ice  conductivity.  When  open  water  absorbs  heat,  the  heat  mixes  underneath  the 
floes  to  decrease  the  vertical  growth  rate.  Any  remaining  heat  either  causes  lateral  melting 
or  raises  the  temperature  of  the  fixed  depth  (30  m)  ocean  mixed  layer.  In  the  presence  of 
an  ice  cover,  the  mixed  layer  temperature  is  always  set  equal  to  freezing.  Thus,  the  excess 
heat  absorbed  by  the  leads  is  used  for  lateral  melting  until  the  ice  disappears.  Also,  during 
growth  conditions,  ice  is  not  allowed  to  form  until  the  mixed  layer  reaches  the  freezing 
temperature  of  sea  water.  For  a  detailed  description  of  the  thermodynamic  portion  of  the 
model  see  Hibler  (1980)  and  Preller  (1985). 

B.  Model  Design 

The  operational  sea  ice  model  is  defined  on  a  subsection  of  the  FNOC  63  x  63  Northern 
Hemisphere  polar  stereographic  grid.  The  model  domain  includes  the  central  Arctic,  the 
Barents  Sea  and  the  Greenland  and  Norwegian  Seas  bounded  at  a  line  extending  from 
Greenland  through  Iceland  to  the  United  Kingdom  (Figure  76).  These  model  boundaries 
are  limited  by  the  available  oceanic  forcing  (see  Hibler  and  Bryan,  1984).  The  model 
dimensions  are  47  x  25  with  127  km  grid  resolution.  All  boundaries  of  the  model  are 
solid  walls  except  for  the  southern  boundary  through  the  Greenland  and  Norwegian  Seas. 
Outflow  grid  cells  have  been  designated  at  this  boundary.  All  ice  which  advects  into  these 
cells  flows  out  of  the  basin. 

The  forcing  for  the  model  may  be  divided  into  two  parts:  atmospheric  and  oceanic. 
Atmospheric  forcing  from  the  NOGAPS  model  consists  of  surface  winds,  surface  air  tern- 


perature,  surface  pressure,  surface  vapor  pressure,  incoming  solar  radiation,  surface  sensible 
heat  flux  and  net  surface  heat  flux.  The  last  three  fields  are  used  to  determine  long  wave 
radiation.  All  of  the  NOGAPS  forcing  fields  are  stereographic  grid  to  the  PIPS  grid.  The 
oceanic  forcing  is  derived  from  the  Hibler/Bryan  ice  ocean  model  (1984).  These  fields  are 
interpolated  from  the  Hibler/Bryan  model  grid  to  the  PIPS  grid. 

C.  Model  Results 

The  model  and  the  forcing  fields  were  tested  via  a  number  of  equilibrium  “spin  up” 
runs.  In  all  of  these  tests  the  model  was  forced  with  the  1983  NOGAPS  fields,  cyclically, 
for  a  three  to  five  year  “spin  up”  period.  The  model  was  initialized  by  setting  the  ice 
drift  equal  to  zero,  defining  the  thickness  field  to  be  3.3  m  everywhere  and  defining  the 
compactness  to  be  1.0  at  all  grid  points.  Table  12  presents  the  parameters  used  in  all  of 
the  test  cases  as  well  as  in  the  final  version  of  the  PIPS  model. 

These  initial  tests  of  the  model  show  that  the  resultant  thickness  and  concentration 
fields  contain  the  observed  seasonal  characteristics  of  ice  concentration  and  thickness.  In 
particular,  the  ice  thickness  fields  show  the  thickest  ice,  approximately  6  -  7  m  (LeSchack. 
1980),  building  up  along  the  Canadian  Archipelago  and  north  Greenland  coast.  Ice  thick¬ 
ness  decreases  to  the  north,  away  from  the  Canadian  Archipelago  towards  the  North  Pole, 
and  continues  to  decrease  from  the  pole  towards  the  Soviet  coast.  These  patterns  exist 
during  all  seasons  with  the  average  ice  thickness  being  smaller  in  the  summer  and  fall  then 
in  the  winter  and  spring  (Figure  77).  Mean  annual  model  derived  ice  thicknesses  (Figure 
78)  were  compared,  by  regions,  to  ice  thickness  derived  from  submarine  data  (Garrett. 
1985).  Results  show  good  agreement  between  model  ice  thickness  and  data.  Although  the 
model  ice  thickness  appears  biased  thin,  the  mean  difference,  0.6  m  is  within  the  error 
of  the  Garrett  data.  Some  difference  may  also  be  due  to  the  fact  that  we  are  comparing 
1983  model  ice  thicknesses  to  three  different  years  of  submarine  data.  It  should  also  be 
noted  that  the  model  thickness  gives  a  good  representation  of  the  relative  variation  in 
observed  ice  thickness  for  the  given  regions.  The  ice  drift  velocities  are  highly  correlated 
to  the  1983  NOGAPS  winds  and  also  agree  well  with  observation.  Another  important 


factor  determined  by  these  tests  is  the  strong  effect  of  oceanic  heat  flux,  particularly  in 
the  Marginal  Ice  Zone.  Tests  which  compared  the  results  using  a  constant  heat  flux  of  2 
watts  per  square  meter,  a  ralid  estimate  for  the  central  Arctic,  to  the  results  using  the 
Hibler/Bryan  monthly  mean  values  showed  a  much  improved  ice  edge  in  the  second  case. 
Oceanic  heating  appears  to  be  a  very  important  factor  for  melting  ice  in  the  Barents  Sea 
and  the  Greenland  Sea. 

The  model  predicts  large  seasonal  variability  in  the  marginal  ice  zone.  A  typical 
winter  thickness  field  shows  the  southern  half  of  the  Barents  Sea  to  be  basically  ice  free 
with  a  maximum  of  approximately  2.5  -  3.0  m  of  ice  piling  up  against  Spitzbergen.  Model 
results  also  show  the  Barents  Sea  to  be  almost  totally  ice  free  in  the  summer.  Along  the 
Greenland  coast,  the  ice  extends  south  to  cover  the  northern  coast  of  Iceland  in  the  winter 
but  retreats  dramatically  in  the  summer. 

D.  Operational  Model 

The  PIPS  model  is  presently  undergoing  its  OPTEST  (validation  by  the  users  of  the 
product).  This  OPTEST  consists  of  two  phases.  In  Phase  I,  the  model  is  initialized  by 
its  own  24  hour  forecast  from  the  previous  day.  If  the  24  hour  forecast  is  not  available, 
the  model  is  initialized  from  climatology.  This  climatology  is  a  model  derived  climatology 
calculated  using  the  1983  NOGAPS  atmospheric  forcing  and  the  Hibler/Bryan  oceanic 
forcing. 

The  first  half  of  the  OPTEST  showed  good  agreement  between  the  predicted  ice  drift 
and  observation  (Tucker  and  Hibler,  1986).  The  model  ice  thickness  and  concentration 
showed  good  agreement  with  the  observed  trends  of  ice  growth  and  retreat.  However,  the 
model,  when  strickly  forced  by  the  NOGAPS  forcing  and  the  Hibler/Bryan  ocean  forcing, 
predicts  too  much  ice  in  the  Barents  Sea  and  the  Greenland  Sea  in  the  winter.  We  feel 
that  these  problems  will  be  alleviated  in  the  coming  year  when  we  include  the  temporal 
variability  of  the  oceanic  forcing  by  including  a  more  sophisticated  ocean  model  under  the 
ice. 
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The  ice  edge  problem  is  presently  being  addressed  in  the  second  phase  of  the  OPTEST. 
In  Phase  II,  the  PIPS  model  is  updated  (initialized)  with  data.  The  data  used  for  updating 
are  ice  concentration  values  derived  from  the  N aval  Polar  Oceanography  Center’s  (NPOC) 
weekly  analysis  of  ice  concentration.  NPOC  creates  a  subjective  hand  analysis  of  ice 
concentration  for  both  the  eastern  and  western  Arctic.  The  analysis  is  derived  from  various 
sources  of  remotely  sensed  data  (infrared  and  passive  microwave)  as  well  as  any  available 
ship  and  aircraft  observations.  The  data  is  digitized,  sent  to  FNOC  and  then  interpolated 
to  the  ice  model  grid.  Once  per  week,  the  model  concentration  field  is  entirely  replaced  by 
the  NPOC  concentration  field  with  the  necessary  adjustments  being  made  to  the  thickness 
fields  and  heat  absorption  in  the  mixed  layer.  Initial  PIPS  model  results  which  have  been 
updated  from  the  NPOC  analysis  show  substantial  improvement  in  the  model’s  ability  to 
predict  ice  edge  location. 

The  PIPS  model  produces  a  144  hour  forecast  (determined  by  the  present  NOGAPS 
forecast  capability)  on  a  daily  basis.  On  days  when  the  model  is  not  being  updated  from 
the  NPOC  analysis,  it  is  initialized  by  the  previous  day’s  24  hour  forecast.  The  operational 
model  outputs  the  following  11  fields  daily: 

1.  Ice  thickness;  Tau  0,  Tau  144 

2.  Ice  concentration;  Tau  0,  Tau  144 

3.  Divergence/convergence;  Tau  48,  Tau  96,  Tau  144 

4.  Ice  displacement  (drift);  Tau  24,  Tau  144 

5.  Six  day  change  in  ice  thickness 

6.  Six  day  change  in  ice  concentration 

where  “Tau”  indicates  the  forecast  time  in  hours.  These  fields  will  be  used  by  the  Naval 
Polar  Oceanography  Center  as  guidance  in  the  generation  of  a  variety  of  products. 
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III.  CONCLUSIONS 
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The  Polar  Ice  Prediction  System  (PIPS),  based  on  the  Hibler  (1979)  ice  model,  has 
the  capability  of  predicting  fields  of  ice  concentration,  thickness  and  drift.  It  has  been 
shown  that  these  three  model  derived  fields  show  good  agreement  with  observations.  As 
an  operational  forecast  tool,  the  PIPS  model  will  predict  these  fields  out  to  144  hours. 
Updating  with  data,  in  the  form  of  digitized  fields  of  ice  concentration  from  the  NPOC 
hand  analysis,  results  in  a  marked  improvement  in  the  model’s  ability  to  predict  ice  edge 
location. 
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Appendix  A 


An  Experiment  in  Data  Assimilation 
Using  the  Kalman  Filter 


Section  _1 .  Introduction 

This  is  an  outline  of  an  experiment  currently  underway  studying  the  use 
of  the  Kalman  Filter  for  assimilating  data  into  a  one  dimensional  model  of 
the  mixed  layer.  We  will  use  the  non-advective  Mellor-Yamada  Level  2  model 
(MYL2)  initialized  from  an  XBT  profile  at  OWS  Romeo  to  make  3  and  5  day 
forecasts.  Approximate  sea  surface  temperatures  from  bucket  thermometers 
will  be  assimilated  into  the  model  forecasts  as  they  become  available.  The 
resulting  statistical-dynamical  forecasts  will  be  compared  to  observed  XBT 
profiles,  and  the  RMSE  errors  compared  to  those  of  persistence  and  of  the 
model  without  data  assimilation. 

There  are  two  purposes.  First,  we  will  study  the  ability  of 
approximate  surface  temperature  information  to  keep  a  model  forecast  'on 
track'  at  a  variety  of  depths  in  the  mixed  layer  and  in  the  seasonal 
thermocline.  Second,  we  will  illustrate  the  many  problems  involved  in  an 
application  of  the  Kalman  Filter  to  a  useful  physical  model,  and  document 
the  ways  those  problems  are  approached  in  this  experiment.  This  should 
provide  a  set  of  experience  for  future  applications. 

I  begin  by  outlining  the  mechanics  of  the  Kalman  Filter,  and 
enumerating  some  of  the  technical  problems  in  its  use.  In  Section  3,  1 
define  the  Filter  for  this  application.  In  Section  4,  I  discuss  proposed 
solutions  to  the  technical  problems. 
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Section  II.  The  Kalman  Filter 

We  begin  by  assuming  the  dynamic  model  can  be  written  as  a  linear 
system. 


k+1 


*k+lwk  +  Fk+luk+l  +  ek+l 


(2.1) 


where  w^  is  an  nxl  vector  denoting  the  true  state  of  the  system  at  time  k. 
Tk+1  is  an  nxn  transition  matrix  expressing  the  dynamics  of  the  system* 
whose  elements  are  independent  of  the  state  vector  w^.  uk+1  is  an  mxl 
vector  of  forcing  conditions,  and  an  nxm  matrix  expressing  their  effect 
on  the  new  state  vector.  We  make  the  simplifying  assumption  that  F  and  u 
are  known  exactly. 

We  assume  that  is  an  nxl  random  vector  with  =  0, 
var(ek+1)  =  Qk+1*  and  cov(£k' ek+j*  =  0  if  )*0.  is  the  error  covariance 
matrix  expressing  the  mean  square  error  in  the  model  when  initialized  from 
perfect  data  and  forecast  forward  a  single  time  step. 

Beginning  with  an  initial  estimate  of  the  state  of  the  system  w£  and  an 
error  covariance  matrix  for  this  estimate,  Pa,  we  proceed  as  follows.  At 
each  time  step: 


T  wa  +  F  u 
Tk+1  k  rk+luk+l 


(2.2) 


wk+l  is  the  dynamic  forecast  for  time  k+1.  It  is  made  by  projecting 
forward  the  analysis  w^  for  time  k.  The  expected  squared  error  for  the 
forecast  is  described  by  the  error  covariance  matrix  p£ 


Pk+1  ~  5  ^  (wk+l'  wk+l>(wk+l'  wk+l}  ^  =  Tk+lPk  Tk+1+  Qk+1  (2,3) 


Here  Pa  is  the  error  covariance  matrix  for  the  analysis  w^. 


c  (  (“k  •  “kl<uk 


(2.4) 


Let  be  a  pxl  vector  of  observations  obtained  since  time  k.  We  assume 
that  these  observations  are  related  to  the  true  state  vector#  and  the 
relationship  is  expressed  by: 


xk+l=  Hk+lwk+l  +  nk+l 


(2.5) 


where  H^+1  is  a  known  pxn  matrix#  and  a  pxl  random  vector  with  mean  0 

and  covariance  matrix  We  define  an  innovations  vector: 


zk+l  xk+l  “  Hk+lwk+l 


(2.6) 


2k+1  expresses  the  apparent  discrepancy  between  the  model  forecast,  w£+1  and 
the  observations.  When  z^+i  is  0,  no  discrepancy  was  observed.  The  Kalinan 
Filter  blends  the  innovations  with  the  forecast  using  the  Kalman  Gain  Matrix 


Pk+lHk+l  ^  Hk+lPk+lHk+l  +  Rk+1  ^ 


f  T 
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The  analysis#  which  is  a  blend  of  model  forecast  and  data  is 


wk+l  wk+l  +  Kk+lzk  =  wk+l+  Kk+l(xk+l“  Hk+lwk  ) 


The  error  covariance  matrix  for  the  analysis  is  computed  as 


(2.7) 


(2.8) 


Pk+1  ^  1  "  Kk+lHk+l)  Pk+1 


(2.9) 


At  the  end  of  this  iteration#  a  new  forecast  is  begun  for  the  next  time 
step  using  wj^  and  Pj^j  which,  together#  describe  the  filter  state  at  any 


one  time. 


Any  application  of  the  Kalman  Filter  first  requires  that  several 
numerical  issues  be  settled/  and  various  covariances  used  by  the  Filter  be 
estimated.  Here  we  list  these  items  where  the  reader  can  observe  their  part 
in  Equations  2. 1-2. 9.  In  Section  4,  we  discuss  the  approach  used  in  this 
work,  for  each  item. 

1.  Actual  dynamic  models  are  rarely  linear.  The  transition  matrix 
often  contains  components  which  are  functions  of  the  state  variables. 

2.  The  bulk  of  the  computational  load  is  in  the  update  of  the  forecast 
covariance  matrix  (Eqn  2.3).  Each  multiplication  of  T  by  a  column  of  P  is 
computationally  equivalent  to  an  iteration  of  the  forecast.  Since  2n  of 
these  multiplications  are  performed  in  each  time  step,  it  is  easily  seen 
that  this  work  must  be  made  as  efficient  as  possible  for  the  computation  to 
be  feasible. 

3.  To  begin  the  iteration,  we  require  an  initial  estimate  of  the  state 
vector,  Wq,  and  an  estimate  of  its  error  covariance  matrix  P^. 

4.  To  compute  the  new  error  covariance  matrix  for  the  forecast,  p£,  we 

need  the  error  covariance  matrix  Q^.  This  matrix  describes  the  error 

properties  of  the  model,  initialized  from  perfectly  correct  data,  for  a 

single  time  step.  While  the  entries  in  this  matrix  are  usually  small,  they 

have  an  important  cumulative  effect  on  the  Filter.  It  is  easily  seen  that 

if  P^=0  and  Qk=0  ,  then  p£=0.  The  fil*er  will  assume  the  model  is 

perfect,  and  discard  the  information  in  the  observations  xk  when  computing 

the  analysis.  As  0^  increases,  the  model  is  assumed  decreasingly  accurate, 

.  » 

and  more  weight  will  be  placed  on  x^  in  the  analysis. 

5.  The  accuracy  of  the  observations  x^,  as  expressed  in  Rk  must  be 
estimated.  The  relative  size  of  R^  and  entries  in  P^  control  the  weight 
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given  the  observations  relative  to  the  model  forecast  w[  in  computing  the 
analysis  w^. 


Section  III .  Experimental  Design,  An  Overview 

The  non-advective  form  of  the  MYL2  model  forecasts  values  of 
temperature  (T)  and  the  Ekman  velocities  (U  and  V)  at  selected  depths  on  a 
numerical  grid.  On  a  grid  with  N+l  points,  the  state  vector  has  dimension 


Tkm 


T.  (N) 


Vk(N) 


(3.1) 


The  values  at  the  N+ls  depth  are  held  constant  throughout  the  forecast 

-« 

period,  so  are  not  included  in  the  state  vector. 

Suppose  an  approximate  sea  surface  temperature  is  available  at  time 


xk+,=  SST  at  time  k+1  =  (1  0  0  ...  0)  f  T. . _ (1)1  +  n. 


Hk+lwk+l  +  nk+l 


(3.2) 


Hk+1  is  the  lx3N  row  vector  with  a  1  in  the  first  place  and  zeroes 
thereafter,  nk+1  is  a  scalar  which  represents  the  difference  in  the  bucket 
temperature  and  the  true  temperature. 

At  OWS  Romeo  (47°N,  17°W)  XBT  temperature  profiles  are  available  every 
12  hours.  Surface  observations  including  sea  surface  temperature'  from 


buckets  and  heat  fluxes  are  available  every  three  hours.  (Some  of  this 
data  is  missing.)  We  initialize  the  model  with  an  XBT  temperature  profile 
on  the  interval  from  0  to  200  meters.  A  typical  short  term  forecast  uses 
time  steps  of  600  seconds  and  a  grid  spacing  of  5  meters.  N  is  40,  and  the 
state  vector  has  dimension  nxl  =  120x1.  Every  3  hours,  or  18  time  steps, 
an  approximate  SST  will  be  blended  with  the  forecast  using  the  Kalman  Filter 
as  described  in  Section  2.  The  RMSE  of  the  three  and  five  day  forecast 
errors  for  temperature  will  be  compared  to  those  of  persistence,  and  to 
those  of  the  model  forecast  without  using  the  assimilation  of  the  SST. 

This  scheme  is  summarized  in  Figure  79. 


In  this  Section/  I  discuss  the  ways  Issues  1-5/  listed  at  the  end  of 


Section  2,  have  been  or  will  be  decided. 

1.  Non-linearity  in  the  model 

The  MYL2  model  is  non-linear  in  that  the  transition  matrix  uses 

entries  which  depend  on  eddy  coefficients/  which  in  turn  depend  on  the 
temperatures  and  velocities.  For  example#  the  updating  procedure  for  T  is 
def ined  by 


Lk+lTk+l  1k  +  Bk+1 


(4.1) 


Dk+1  is  an  nxl  vector  expressing  surface  forcing  and  the  extinction  of 
solar  radiation.  L  is  a  tridiagonal  matrix  whose  coefficients  are  functions 
of  the  Richardson  number  computed  using  T^#  and  V^. 

The  usual  approach  for  nonlinear  models  is  the  Extended  Kalman  Filter. 
In  this  method#  the  entries  in  -L  (or  Y)  are  simply  computed  as  estimated 
from  the  current  analysis  of  the  state#  and  the  resulting  errors  in  the 
computation  of  the  error  covariance  matrix  are  ignored.  This  is  the  method 
adopted  here. 

We  will  ignore  errors  in  the  components  of  L»  essentially  treating  them 
as  if  they  were  known  perfectly.  This  necessarily  introduces  some  sub-^ 
optimal  performance  into  the  filter,  and  some  underestimation  of  the  true 
variance  of  the  forecast  error.  The  effect  of  this  and  other  filter  design 
imperfections  will  be  seen  when  we  compare  the  size  of  actual  errors  in  the 
statistical  dynamical  forecasts  with  the  computed  expected  error 
•covariances. 


2.  Computation  of  p£  explicitly  or  implicitly. 


As  written  in  Eqn  2.1,  the  Kalman  Filter  requires  an  explicit 


representation  using  a  transition  matrix  T^.  In  fact,  the  numerical  code 
used  at  NORDA  (written  by  P.  Martin)  uses  an  implicit  scheme  as  written  in 
Eqn  4.1.  There  are  two  choices'  for  procedure. 

a.  An  explicit  representation  can  be  achieved  from  (4.1)  by  inverting 
the  tridiagonal  matrix  L^+j.  This  produces 


Tk+i  =  LkliTk  ♦ 


Lk+1  °k+l 


=  JL.,T.  + 

k+1  k  k+1  k+1 


(4.2) 


The  error  covariance  p£+j,  produced  as  in  Eqn  2.3  requires  two 
multiplications  of  nxn  matrices  to  produce  Jr+jPrTr+i  •  This  time 
consuming  process  is  best  avoided. 

b.  The  new  covariance  matrix  p£+1  can  be  produced  using  the  implicit 

representation  with  If  P^)  is  the  i^  column  of  P^,  then  a 

multiplication  P*^  =  *R+ip(i)  equivalent  to  solving 

Lk+lP(i>  =  P(i)‘  Hence,  the  product  Jr+iPrTr+i  may  be  produced  by 

first  applying  the  implicit  algorithm  n  times  to  columns  of  P  to 

produce  1r+ipr*  A  transpose  produces  Pr1r+i'  and  another  n 

a  T 

applications  of  the  implicit  procedure  produces  *k+lPkTk+l  * 

The  implicit  procedure  described  in  b  is  considerably  faster  than  the 
explicit  multiplication.  However,  the  explicit  procedure  need  not  be  fully 
carried  out  at  each  time  step.  Instead,  the  cumulative  transition  matrix 
A  =  Tr+itr1r-i*  •  •  need  only  be  updated  at  each  step,  and  the  matrix  p£+1  be 
produced  only  when  observations  are  to  be  blended. 

Both  an  implicit  and  explicit  scheme  have  been  coded  and  checked  at  the 
University  of  North  Florida  by  D.  Mohr.  This  code  will  be  shipped  to  J. 
French  at  NORDA  by  Oct.  5,  1986.  After  it  is  installed  and  rechecked  at 
NORDA,  benchmark  timing  tests  will  be  run  to  compare  the  speed  of  each 


algorithm. 

3.  Initialization  of  w®  and  P^. 

The  forecasts  are  initialized  by  setting  the  temperatures  equal  to 
those  of  an  interpolated  XBT  profile  and  the  velocities  equal  to  0.  He  will 
assume  there  is  no  error  in  the  temperature  Initialization. 

The  error  symmetric  covariance  matrices  P  can  be  better  understood  by 
looking  at  the  following  decomposition: 

pTT  pTU  pTV 

P*  =  pUT  pUU  pUV 
pVT  pVU  pW 

Here  P^  is  the  covariance  matrix  for  the  errors  in  temperature#  P^U  is 
the  cross  covariance  matrix  for  the  errors  in  temperature  with  those  in 
velocity#  etc.  The  assumption  that  there  is  no  error  in  the  initialization 
for  temperature  implies  pTri=pTU=pTV_pUl_pVT_o 

To  estimate  reasonable  values  for  the  error  variance  in  the  velocities, 
we  turn  to  the  MILE  Data.  This  data  set  contains  temperatures  and 
velocities  for  a  series  of  observations  at  a  station  PAPA  in  the  North 
Pacific.  While  this  Station  is  very  different  from  OHS  ROMEO#  its  data  can 
be  used  to  at  least  estimate  the  magnitude  and  shape  of  covariance 
functions. 

The  entries  in  PUI^#  for  example#  are  &  (U(i)U()))  where  U(i)  is  the 
actual  U  velocity  at  depth  i.  He  turn  to  the  sample  moments  observed  in  a 
subsample  of  134  observations  in  the  MILE  data. 

uvi5  =  J  U(i)V(j)  /  134 

The  sample  moments  uu^  and  vv^  have  maximum  values  of  27#  000  and 


19/000  cm  /sec  at  the  surface.  The  moments  decrease  to  around  10.000  and 
7/000  cm^/sec^  by  the  50m  depth.  As  a  first  approximation  for  the  entries 
along  the  diagonals  of  P^  and  P^  1  propose  a  value  of  25.000  in  the 
initial  position  (for  depth  OiiO,  decreasing  linearly  to  9. 000  at  the  11^ 
position  (50m)  and  remaining  constant  thereafter. 

Errors  at  nearby  depths  are  surely  closely  related.  However, 
initially,  we  will  take  values  off  the  diagonals  in  Puu  and  Pw  to  be  zero, 
together  with  all  the  entries  in  P1^  and  P^°.  By  examining  the  cross  sample 
moments  in  the  MILE  Data,  these  entries  will  be  better  estimated  in  the 
future. 

4.  Estimation  of  Q^. 

To  estimate  the  error  properties  of  the  model,  initialized  from 
•perfect'  data,  we  turn  to  the  MILE  Data  mentioned  in  Part  3  above.  A 
subsample  of  the  large  data  set  was  created  by  selecting  pairs  of  records 
one  hour  apart  in  time,  with  the  first  record  of  each  pair  7  hours  apart  in 
time.  That  is,  records  were  selected  at  times  0,1,7,8,14,15,...  hours.  The 
goal  of  selecting  pairs  7  hours  apart  was  first  to  obtain  observations  more 
nearly  independent  than  those  112.5  seconds  apart  as  in  the  original  data 
set.  Second,  by  selecting  a  7  hour  stagger,  our  sampling  interval  is 
relatively  prime  to  the  24  hour  diurnal,  12  hour  semidiurnal  and  15-16  hour 
inertial  frequencies. 

For  each  pair,  the  model  was  initialized  with  the  temperatures  and 
salinities  for  the  first  record,  and  forecast  foreward  for  six  600-second 
time  steps,  or  one  hour.  The  resulting  forecast  was  compared  to  the  actual 
observations  in  the  second  forecast.  The  sample  covariances  for  the  errors 
have  been  computed  by  Jon  French.  As  of  9/25/86  we  have  begun  smoothing 
this  rough  data,  to  find  a  tractable  approximate  covariance  function. 


IJI  JUVIWlA’PMfl'J 


Riivivv  vvvtvi  n  rwvm 


■t 


E 


f 


S 


I  r 

•  * 

« 

*  * 
r.- 

i 


5.  Estimation  of  R^. 

Since  our  observations  consist  of  a  single  number,  the  sea  surface 
temperature  from  bucket  readings,  R^  consists  of  simply  the  Mean  Square 
difference  in  the  bucket  temperature  and  the  temperature  from  an  XBT  near 
the  surface,  within  one  time  step  (600  seconds)  of  each  other.  From 
historical  data  at  OWS  Romeo,  we  know  that  a  sea  surface  temperature  from 
buckets,  and  the  shallowest  point  on  an  XBT  thermal  profile  taken  an  hour 
later,  will  have  a  Root  Mean  Square  Difference  of  about  .16°C.  This  gives  a 
maximum  value  of  R^  of  (.16)2.  Since  the  Filter  will  use  time  steps  of  10 
minutes,  we  assume  an  appropriate  value  for  R^  is  more  like  (.05)2  to 
(.10)2. 


Appendix  B 


1  .  Statistical  Methods 

A  process  Z(t)rx^  consisdng  of  r  measurements  at  time  t  is  said  to 
follow  a  stationary  mean  0  vector  autoregressive  process  of  order  1  If 


E(Z(t))  -  tfrxl 

(1) 

G(k)rxr  “  cov(Z(t+k),Z(t)) 

for  all  t 

(2) 

Z(t+1)  -  Brxr  Z(t)  + 

e(t) 

(3) 

The  Innovation  process  e(t)  Is  assumed  to  follow  a  multivariate  normal 
distribution  with  mean  0,  variance  matrix  Vrxr  and  cov(e(t+k),e(t))  *  0  when 


Under  these  conditions 

G(k)  -  Bk  G(0) 

k  >  J6 

(4) 

G(0)  -  B  G(0) 

B*  +  V 

(5) 

Both  G(0)  and  V  must  be  positive  definite  matrices,  if  the  statlonarity 
condition  In  (2)  is  fulfilled. 

In  an  observed  series  0(t)  where  persistence  Is  very  strong,  the 
statlonarity  condition  often  fails.  Then  the  observed  process  Is  usually 
differenced  so  that  Z(t) *  0(t)  -  O(t-l). 

Once  B  and  G (fi)  (or  V)  are  estimated,  linear  statistical  forecasting 
proceeds  via  the  Gauss-Markov  Theorem  (see  Bretherton,  Fandry  and  Davis  1976.) 
Let  Y  be  a  vector  to  be  predicted,  and  X  a  vector  of  predictor  variables, 
both  with  mean  0.  Then  the  unbiased  minimum  variance  linear  prediction  of  Y 
Is  Y'  -  X'  E(XX'rlE(XY'). 

In  our  application  0'(t)  -  (0j'(t)l02'(t»  consists  of  Oj(t),  the  vector 
of  3  EOF  coefficients  used  to  represent  deviations  from  climatology  at  time  t, 
and  02(t)  is  a  vector  of  measurements  describing  surface  conditions  at  time  t. 
Z(t)  “  0(t)  -  O(t-l)  Is  the  vector  representing  changes  In  deviations  and 


changes  In  surface  conditions 


We  wish  to  predict  Y  ■  Oj(t+k)  -  Oj(t)  ■  ZZ|(t+i)  given 

zx(t) 

OjCt-HO-OjCt) 

The  covariance  matrices  E(XX')  and  E(XY')  can  be  computed  from  B  and  G(0) 
using  (4).  The  predicted  temperature  profile  can  be  constructed  from  Y  by 

Profile(t+k)9xl  -  Eg  Oj(t)  +  Eg  Y  +  climatology(t)<)xl 

where  Eg  is  the  9x3  matrix  whose  columns  are  the  EOF's  for  appropriate  season 
s. 

2 .  Estimating  B,  G(0) 

First  estimates  of  B  and  G(j3),  denoted  B*  and  G*  are  provided  by  the 
ordinary  multivariate  regression  of  the  observed  Z(t)  on  Z(t-l)  and  the  lag 
(0)  sample  covariance  matrix  of  the  Z(t).  Two  problems  arise.  First,  because 
ocean  data  sets  are  incomplete,  many  observations  will  not  be  paired  with  an 
observation  on  the  previous  day.  A  method  which  only  uses  points  where 
immediately  preceding  observations  exist  wastes  valuable  data.  Second,  the 
resulting  B*  and  G*  will  not  necessarily  represent  a  stationary  process,  since 
V*  *  G*  -  B*G*B*#  (applying  (5))  will  not  necessarily  be  positive  definite. 
The  fitted  process  will  be  useless  for  Inference  or  simulation. 

The  maximum  likelihood  criterion  will  be  applied  to  the  original 
estimates.  The  negative  log-likelihood  of  the  observed  Z(t)  at  times'^ 

t2*t2» . tn  conditional  on  Z(tj)  is  approximately  minimized  using 

Marquardt-Levenberg  procedures,  first  over  choices  of  B,  leaving  G  fixed  at 
G*,  then  over  choices  of  G,  and  finally,  over  B  again.  By  forcing  the 
minimizing  routine  to  compute  a  large  value  when  the  choice  of  B  or  G  yields 


'  ZjCt) 

Z  Z~(t+i) 


non-positive  definite  V,  the  fitted  process  will  be  an  acceptable  one. 

3 .  Avoiding  Overfitting 

The  dangers  of  overfitting  in  oceanographic  data  have  been  discussed  by 
Davis  (1978),  Barnett  and  Hasselmann  (1979)  and  others.  We  shall  control  for 
overfitting  by  cross-validation,  that  is  by  checking  the  results  developed 
using  one  year's  season  on  the  data  from  the  other  year's  season.  This  method 
has  been  found  superior  to  all  others,  since  it  allows  the  researcher  to  check 
for  shifts  in  the  underlying  physical  process  which  would  nullify  any  apriori 
computation  of  skill.  (See  Berk,  1984.) 

When  separate  seasons  are  not  available  for  cross-validations, 
researchers  require  some  calculation  of  skill  which  will  guide  them  in  the 
selection  of  informative  covariates  (our  surface  information,  or  02(t).)  In 
the  autoregressive  setting  in  which  we  are  working,  the  most  appropriate  is 
the  Akalke  Final  Prediction  Error  (Akaike,  1971.)  Let  S(M)  be  the  sample 
covariance  matrix  of  the  lag-one  prediction  errors  from  a  model  M  (that  is,  a 
list  of  covarlates  02(t)  and  estimates  of  B  and  G(0).  Akaike  recommends 
choosing  M  to  minimize 

MFPE(M)  -  (l-(P4tj)/n)"P  (l+(p+q)/n)P  S(M)  (6) 

Here  p  is  the  dimension  of  the  process  being  predicted  (in  our  case  3)  and 
q  the  number  of  covarlates,  that  is  the  dimension  of  02(t).  The  determinant 
S(M)  will  be  large  when  models  are  underfit,  and  will  reduce  as  covarlates  are 
added.  The  first  two  factors  will  Increase  as  q  «  number  of  covarlates  is 
Increased.  The  minimum  is  achieved  when  goodness  of  fit  is  balanced  off 
against  the  penalty  for  extra  covariates.  In  practice,  Akaike's  method  yields 
a  few  models  with  similar  MFPE  all  of  which  should  give  acceptable  results. 
We  will  check  the  results  from  the  HFPE  analysis  against  the  results  from  the 
cross-validations.  Large  discrepancies  in  the  results  must  mean  that 
forecasts  developed  in  one  year  are  not  always  appropriate  the  next. 
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FIGURE  1:  SHALLOW  OCEAN  GRID 
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FIGURE  2:  DEEP  OCEAN  GRID 
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Figure  3.  Plane  wave  incident  from  fluid  1  to  fluid  2,  with  c0<d  and  a  >  <f>- 
c0  is  sound  speed  in  fluid  1,  c7  is  sound  speed  in  fluid  2,  a  is  angle  of  incidence 
(reflection)  and  <j>  is  angle  of  transmission. 
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Figure  4.  Simulation  window  with  initial  waveform.  Regions  1  and  2  are  fluids  with  sound 
speeds  1500  m/s  and  1600  m/s  and  depths  of  400  m  and  80  m  respectively.  Region  3  is  a 
damping  region  with  sound  speed  1600  m/s  and  depth  of  120  m.  The  width  is  200  m. 
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Figure  7.  Simulation  window  with  initial  waveform  for  deep  ocean.  Initial  range  of  frame 
is  4194  m  (range  of  left  hand  side).  Width  of  frame  is  1050  m,  total  depth  of  frame  is 
6000  m  with  bottom  400  m  for  damping.  Sound  speeds  in  ocean  are  1540  m/s  at  surface, 
1495  m/s  at  depth  of  1000  m  and  1560  m/s  at  depths  of  5500  m  to  6000  m.  Intermediate 
depths  are  linearly  interpolated  from  these  values. 
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Figure  8.  Linear  deep  ocean  simulation  (fi  =  0)  with  initial  conditions  of  Figure  7. 


RANGE  - 


•4594.00  M  NSTEP  -  8000 


RANGE  -  78144.00  M  NSTEP  -  MOO 


Figure  9.  Nonlinear  deep  ocean  simulation  (0  =  3.5)  with  initial  conditions  of  Figure  7. 
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Figure  10.  Linear  ice  simulation  (/3  =  0). 
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Figure  13.  Linear  case.  Hydrophone  readings  at  (a)  600  m  and  (b)  20000  m  in  range,  and  (c)  mode 
structures  found  at  range  of  20000  m.  Dashed,  vertical  lines  in  (b)  mark  where  mode  structures  in 
were  taken. 
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Figure  14.  Nonlinear  case  of  Figure  13. 
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Figure  15.  Linear  case.  Spectra  at  ranges  of  (a)  600  m,  (b)  5000  m,  (c)  10000  m  and  (d)  20000  m. 
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Figure  17.  NPE  (solid)  and  analytic  (dashed)  crc»s-sections  of  pressures  on  gn.t 
ranges  of  (a)  .1  km,  (b)  .9  km,  (c)  2.2  km  and  (d)  4.2  km 
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Figure  18.  Error  plot  for  case  shown  in  Figure  17. 
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Figure  19.  NPE  (solid)  and  analytic  (dashed)  cross-sections  of  pressures  on  grid  at 
ranges  of  (a)  7.0  km,  (b)  9.7  km,  (c)  13.1  km  and  (d)  16.7  km. 
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Observed  (soid  Sne)  and  predated  (dashed  ine)  shear  durirtg 
MLE  between  (a)  5  and  8  m  depth,  (b)  8  and  1 1  m  depth,  (c)  1 1 
and  20  m  depth,  and  (d)  20  and  23  m  depth  The  observed  time 
series  of  shear  have  been  fltered  to  eiminate  frequencies  above 
2  cph 
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Figure 
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Observed  (sofd  Sne)  and  predated  (dashed  Sne)  shear  ckaring 
MLE  between  (e)  23  and  26  m  depth,  (f)  26  and  29  m  depth, 
(g)  29  and  32  m  depth,  and  (h)  32  and  35  m  depth 
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Figure  30 


Variance-conserving  power  spectra  of  E-W  component  of  shear  in  the 
seasonal  thermodne.  The  most  prominent  peak  corresponds  to  the 
inertia]  frequency  (.064  cydes/hr)  and  the  peak  just  to  its  right 
corresponds  to  the  semi-dumal  internal  tide  (.0806  cydes/hr). 

The  other  two  prominent  peaks  occur  at  the  location  of  the  first 
harmonics  of  the  inertia]  and  semi-ciumal  tidal  frequencies. 


Figure  32  Monthly- averaged  vertical  temperature  profiles  at  the  4  North 
Pacific  stations  for  the  month  of  February. 


Figure  33  Monthly-averaged  vertical  temperature  profiles  at  the  4  North 
Pacific  stations  for  the  month  of  March. 


gure  34  Monthly-averaged  vertical  temperature  profiles  at  the  4  North 
Pacific  stations  for  the  month  of  April. 
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averaged  vertical  profiles  of  the  temperature  variance 
North  Pacific  stations  for  the  month  of  February. 


Figure  3 


Figure  40  Monthly- averaged  vertical  profiles  of  the  temperature  variance 
at  the  4  North  Pacific  stations  for  the  month  of  April. 


averaged  vertical  profiles  of  the  temperature  variance 
North  Atlantic  stations  for  the  month  of  March. 


if  the  mixed  layer  depth  for  the  months  of 
i  July  at  the  4  North  Pacific  stations  according 
C  criterion. 


Figure  45  Time  evolution  of  the  mixed  layer  depth  for  the  months  of 

February  through  July  at  the  4  North  Atlantic  stations  according 
to  a  at  «  .2*  C  criterion. 


Time  evolution  of  the  mixed  layer  depth  for  the  months  of  February 
through  July  at  the  4  North  Atlantic  stations  according  to  a 
critical  Richardson  number  R1  «  .23  criterion. 


Distribution  of  the  mixed  layer  depth  for  the  months  of  February 
through  April  at  the  4  North  Pacific  stations  according  to  a 
A  T  »  .2°  C  criterion. 


Distribution  of  the  mixed  layer  depth  for  the  months  of  February 
through  April  at  the  4  North  Pacific  stations  according  to  a 
critical  Richardson  number  R1*.23  criterion. 
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Monthly-averaged  vertical  profiles  of  the  Brunt-Valsala 
frequency  (radians  s-1  x  10°)  at  the  4  North  Pacific  st 
for  the  month  of  April. 


56  Monthly-averaged  vertical  profiles  of  the  Brunt-Yalsala 

frequency  (radians  s_1  x  106)  at  the  4  North  Atlantic  stations 
for  the  month  of  March. 


North  Atlantic  stations 


104)  for  the  months  of 
Atlantic  stations. 


Distribution  of  the  25  m  shear  (s_1  x  104)  for  the  months  of 
February  through  April  at  the  4  North  Pacific  stations. 
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x  104)  for  the  months  of 
Atlantic  stations. 


Figure  64  Distribution 


Distribution  of  the  shear  (s'1  x  1(f)  summed  over  all  17  depth 
levels  for  the  months  of  February  through  April  at  the  4  North 
Atlantic  stations. 
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Distribution  of  the  10  m  Richardson  number  for  the  months  of 
February  through  April  at  the  4  North  Atlantic  stations. 


Figure  67  Distribution  of  the  10  m  Richardson  number  for  the  months  of 
February  through  April  at  the  4  North  Pacific  stations. 


Distribution  of  the  21.3  m  Richardson  number  for  the  months  of 
February  through  April  at  the  4  North  Pacific  stations. 


Figure  69  Distribution  of  the  21.3  m  Richardson  number  for  the  months  of 
February  through  April  at  the  4  North  Atlantic  stations. 


Distribution 


Distribution  of  the  Richardson  number  summed  over  all  17  depth 
levels  for  the  months  of  February  through  April  at  the  4  North 
Atlantic  stations. 


Distribution  of  the  5  m  eddy  coefficient  of  momentum  for  the 
months  of  February  through  April  at  the  4  North  Pacific 
stations. 


Distribution  of  the  5  m  eddy  coefficient  of  momentum  for  the 
months  of  February  through  April  at  the  4  North  Atlantic 
stations. 


Figure  7  4  Distribution  of  the  10  in  eddy  coefficient  of  momentum  for  the 
months  of  February  through  April  at  the  4  North  Pacific 
stations. 
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Seasonal  contours  of  Ice  thickness  from  the  PIPS  model  forced  with  1983 
NOGAPS  Atmospheric  forcing  and  Hlbler/Bryan  ocean  forcing.  Contours  In 
1  m  Intervals. 
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Figure  7 8  Mean  annual  ice  thickness  from  the  model  (hatched 
bar)  and  observations  (open  bar)  for  different  regions. 
Numbers  below  each  bar  indicate  the  ratio,  expressed  as 
percent,  of  modeled  to  observed  thickness. 


TABLE  1.  ROOT  MEAN  SQUARE  ERROR  IN  ONE-DAY  PREDICTIONS. 
All  errors  are  computed  on  the  same  days,  which  had  wind  and 
sea  surface  data  available  for  the  statistical  forecasts. 


SPRING 


Depth 

Climatology 

Persistence 

Stat(SS) 

Stat(WS) 

Stat(NC) 

2.5 

.64 

.23 

.20 

.25 

.29 

7.5 

.61 

.21 

.18 

.22 

.27 

12.5 

.60 

.20 

.18 

.21 

.25 

17.5 

.59 

.35 

.20 

.23 

.26 

25.0 

.62 

.35 

.34 

.34 

.31 

32.5 

.76 

.33 

.36 

.36 

.31 

40.0 

.77 

.35 

.34 

.35 

.30 

50.0 

.74 

.25 

.25 

.25 

.23 

62.5 

.62 

.22 

.27 

.26 

.23 

SUMMER 

2.5 

.91 

.41 

.36 

.43 

.46 

7.5 

.90 

.40 

.35 

.41 

.44 

12.5 

.88 

.44 

.37 

.40 

.41 

17.5 

.94 

.56 

.57 

.52 

.52 

25.0 

1.06 

.80 

.79 

.74 

.77 

32.5 

1.31 

1.08 

.93 

.94 

1.03 

40.0 

1.32 

.89 

.83 

.86 

.92 

50.0 

1.23 

.92 

.84 

.87 

.81 

62.5 

1.15 

.61 

.72 

.75 

.70 

RMSE  In  °C 


TABLE  2.  SURFACE  DATA  AVAILABLE  FOR 
USE  IN  PREDICTIONS 


SST  —  Sea  Surface  Temperature  from  bucket  readings 
average  of  observations  at  18  and  21GMT 

WS  —  Wind  Stress 

average  of  observations  at  12,  15,  18  and  21GHT 

S0L2  -  Solar  Radiation 

average  of  observations  at  12  and  15GMT 

LAT2  -  Latent  Flux 

average  of  observations  at  12  and  15GMT 
LAT3  -  Latent  Flux 

average  of  observations  at  18  and  21GMT 

SEN2  -  Sensible  Flux 

average  of  observations  at  12  and  15GMT 

SEN3  -  Sensible  Flux 

average  of  observations  at  18  and  21GMT 

AIR  —  Air  Temperature 

average  of  observations  at  18  and  21  GMT 


TABLE  3.  CORRELATIONS  OF  CHANGES  IN  TEMPERATURES  AT  VARIOUS 
DEPTHS  WITH  THE  SIMULTANEOUS  CHANGES  IN  SURFACE  VARIABLES. 


SPRINGS 


Depth 

SST 

WS 

2.5 

.884 

-.489 

12.5 

.757 

-.414 

25.0 

-.067 

-.110 

50.0 

-.198 

.184 

SOL  2 

LAT2 

LAT3 

.172 

.164 

.088 

.238 

.186 

.079 

.075 

-.050 

.018 

-.077 

.020 

-.074 

SEN2 

SEN3 

AIR 

.137 

.007 

.231 

.194 

-.060 

.042 

.150 

-.098 

-.308 

-.065 

-.056 

-.054 

SUMMERS 


TABLE  4.  ROOT  MEAN  SQUARE  ERRORS  IN  PREDICTIONS  FOR  SPRING  PROFILES, 
AT  DIFFERENT  FORECAST  LAGS  AND  DEPTH  GROUPS. 


Climatology 

Persistence 

Lag 

n 

S 

M 

D 

S 

M 

D 

1 

48 

.61 

.72 

.68 

.21 

.34 

.24 

2 

42 

.60 

.73 

.68 

.35 

.37 

.27 

3 

40 

.57 

.71 

.62 

.42 

.33 

.36 

4 

41 

.53 

.72 

.62 

.53 

.47 

.37 

5 

39 

.55 

.75 

.63 

.57 

.43 

.41 

6 

37 

.54 

.74 

.60 

.63 

.49 

.43 

7 

39 

.59 

.78 

.65 

.67 

.51 

.46 

Stat 

(  SS 

) 

Stat  (  US 

) 

Stat 

C  NC 

:  ) 

Lag 

n 

S 

M 

D 

S 

M 

D 

S 

M 

D 

1 

48 

.19 

.35 

.26 

.23 

.35 

.25 

.27 

.31 

.23 

2 

42 

.22 

.34 

.27 

.36 

.33 

.27 

.35 

.32 

.26 

3 

40 

.25 

.36 

.36 

.42 

.36 

.35 

.42 

.42 

.42 

4 

41 

.35 

.51 

.38 

.56 

.45 

.36 

.55 

.44 

.27 

5 

39 

.39 

.45 

.39 

.58 

.39 

.38 

.58 

.41 

.40 

6 

37 

.45 

.49 

.41 

.63 

.46 

.40 

.64 

.47 

.43 

7 

39 

.57 

.56 

.43 

.71 

.49 

.43 

.69 

.50 

.45 

RMSE  In  °C 


SQRT 

SQRT 

SQRT 


(  average  squared  error  )  at 
(  average  squared  error  )  at 
(  average  squared  error  )  at 


depths  2.5,  7.5,  12.5  and  17.5  m 
depths  25.,  32.5  and  40  m 
depths  50.  and  62.5  m 


TABLE  5.  ROOT  MEAN  SQUARE  ERROR  IN  PREDICTIONS  FOR  SUMMER  PROFILES 
AT  DIFFERENT  FORECAST  LAGS  AND  DEPTH  GROUPS 


k 

Climatology 

Persistence 

!> 

•  ft 

Lag 

n 

S 

M 

D 

S 

M 

D 

ft 

1 

65 

.91 

1.23 

1.19 

.46 

.93 

.78 

2 

61 

.86 

1.24 

1.07 

.50 

1.07 

.93 

1  m 

3 

62 

.86 

1.34 

1.21 

.53 

1.22 

.89 

s  > 

4 

66 

.87 

1.34 

1.16 

.70 

1.27 

.93 

5 

65 

.94 

1.14 

.94 

.84 

.88 

.91 

\  ^ 

6 

63 

.94 

1.09 

.99 

.93 

1.02 

1.08 

V  s' 

7 

68 

.90 

1.27 

1.13 

.91 

1.41 

1.17 

i 

Stat  (  SS  ) 

Stat  (  US 

) 

1 

Stat  (  NC  ) 

ft 

Lag 

n 

S 

M 

D 

S 

M 

D 

S 

M 

D 

ft 

1 

65 

.42 

.85 

.78 

.44 

.85 

.81 

.46 

.91 

.76 

ft 

2 

61 

.50 

1.16 

.96 

.53 

1.10 

.91 

.53 

1.06 

.87  j 

Hi  F 

3 

62 

.51 

1.31 

1.02 

.59 

1.20 

.85 

.59 

1.24 

.86  1 

ft! 

4 

66 

.67 

1.42 

1.18 

.73 

1.28 

.92 

.73 

1.15 

.91 

Si" 

5 

65 

.73 

1.10 

1.22 

.83 

.91 

.90 

.84 

.87 

.89 

V  > 

6 

63 

.78 

1.26 

1.40 

.97 

1.01 

1.07 

.93 

.98 

1.00 

V  .-S 
s* 

7 

68 

.80 

1.62 

1.56 

.88 

1.38 

1.15 

.90 

1.36 

1.15 

RMSE  in  °C 


S,  M  and  D  as  In  Table  4. 


TABLE  6.  ROOT  MEAN  SQUARE  ERROR  IN  PREDICTIONS  FOR  SUMMER  PROFILES 
EXCLUDING  EDDY  IN  1984,  AT  DIFFERENT  FORECAST  LAGS  AND  DEPTH  GROUPS 


Climatology  Persistence 


Lag 

n 

S 

M 

D 

S 

M 

D 

1 

41 

1.00 

1.01 

.90 

.41 

.63 

.56 

2 

40 

.89 

.99 

.68 

.41 

.81 

.65 

3 

40 

.91 

1.01 

.81 

.45 

.87 

.58 

4 

43 

.92 

1.04 

.83 

.71 

.97 

.60 

5 

44 

1.00 

1.03 

.73 

.88 

.97 

.70 

6 

42 

1.02 

.98 

.68 

.99 

1.10 

.81 

7 

44 

.94 

.93 

.71 

.91 

1.14 

.69 

Stat  (  SS 

) 

Stat 

(  WS 

) 

Stat  (  NC  ) 

Lag 

n 

S 

M 

D 

S 

M 

D 

S 

M 

D 

1 

41 

.33 

.58 

.63 

.41 

.64 

.59 

.42 

.59 

.58 

2 

40 

.35 

.73 

.66 

.46 

.79 

.69 

.46 

.77 

.65 

3 

40 

.40 

.87 

.63 

.51 

.87 

.58 

.51 

.87 

.54 

4 

43 

.57 

1.02 

.82 

.74 

.96 

.68 

.73 

.96 

.64 

5 

44 

.70 

1.15 

.88 

.88 

.98 

.70 

.89 

.97 

.66 

6 

42 

.78 

1.35 

.99 

.99  1 

.08 

.81 

1.00 

1.07 

.72 

7 

44 

.71 

1.37 

.80 

.88  1 

.19 

.77 

.91 

1.19 

.64 

RMSE  In  °C. 


S,  M  and  D  as  In  Table  4 
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TABLE  7.  AKAIKE'S  MFPE  FOR  MODELS  FIT  TO  INDIVIDUAL  SEASONS 


Spring  1983 

Spring  1984 

Stat 

(SS) 

.0003 

.000007 

Stat 

(WS) 

.0113 

.000495 

Stat 

(NC) 

.0075 

.000516 

Summer  1983 

Summer  1984 

Stat 

(SS) 

.025966 

. 164246 

Stat 

(WS) 

.166347 

.414407 

Stat 

(NC) 

.44303 

.72612 

I 


o 


•!:ii 


TABLE  8 A: 

PACIFIC  STATIONS 

Lat/Lon 

t2 

t4 

Bering  Sea 

50N/165E 

- 

5/23 

Papa 

50N/145W 

- 

5/18 

November 

30N/140W 

- 

5/15 

Victor 

34N/164E 

4/7 

5/10 

Mean 


TABLE  B:  ATLANTIC  STATIONS 


Lat/Lon  tl 


Lima 

57N/20W 

3/15 

Charlie 

53N/36W 

3/25 

Romeo 

47N/17W 

3/28 

4/22  5/5 


6/13 


TABLE  9«:  PACIFIC  STATIONS 


TABLE  1  la:  PACIFIC  STATIONS 


Table  1  2  Parameters  used  in  the  numerical  simulation 


